«D-A0<»3  733 


unclassified 

I OF  I 

^043733 


FRANK  J SEILER  RESEARCH  LAB  UNITED  STATES  AIR  FORCE  —ETC  F/6  12/1 
A STUDY  OF  tunable  INTEGRATION  AND  CONTROL  THEORY  FOR  THE  ANALY— ETC(U) 
may  77  ml  SABIN 

FJSRL-TR-77-0006  m 


MICROCOPY  RlkllUTlON  U bl  C HM^l 

NA-i:  HURO*ii  .'•••%  ■ ••  * 


J 


FRANK  J.  SEILER  RESEARCH  LABORATORY 


SRL-TR-77-0006 


MAY  1977 


A STUDY  OF  TUNABLE  INTEGRATION 
AND  CONTROL  THEORY  FOR  THE  ANALYSIS 
OF  DIFFERENTIAL  EQUATION  SOLVERS 


FINAL  REPORT 


PROJECT  230A 


APPROVED  FOR  PUBLIC  RELEASE; 
DISTRIBUTION  UNLIMITED. 


AIR  FORCE  SYSTEMS  COMMAND 


UNITED  STATES  AIR  FORCE 


UNCLASSIFIED 


security  classification  of  This  page  (When  Data  Bntarcdj  

REPORT  DOCUMENTATION  PAGE  , before  completing  form 

"i  REPORT  number  ^ ’,2  GOVT  ACCESSION  NO.l  3.  RECIPIENT’S  CATALOG  NUMBER 


pi  REPORT  number 


rL‘^SRL-TR-77-0006  AD  A- 

4.  TITLE  fand  Srjb(i(/e)  ' 


A STUDY  OF  TUNABLE  INTEGRATION  AND  CONTROL 
THEORY  FOR  THE  ANALYSIS  OF  DIFFERENTIAL 
EQUATION  SOLVERS 


17  AUTHORri) 


Marc  L. 


.^Sabin 


9 PERFORMING  organization  NAME  AND  ADDRESS 

Frank  J.  Seiler  Research  Laboratory  (AFSC) 
USAF  Academy,  Colorado  80840 

n.  controlling  optice  name  and  address 


/ Final  Repiirt,^' 

Jantwir  75  — ‘ DecMhrr  76 


6.  PERFORMING  ORG.  REPORT  NUMBEft 


8.  CONTRACT  OP  GRANT  NUM8ERr«; 


10.  PROGRAM  element,  project.  TASK 
AREA  & WORK  UN  T NUMBERS 


0f^6  61102F 
2304-F1-55 

12.  REPORT  DATE 


Frank  J.  Seiler  Research  Laboratory  (AFSC)  IViL-NT^VR Spaces 
USAF  Academy,  Colorado  80840  7, 


14  MON 'TOPING  AGENCY  NAME  d ADORESSr//  ditterent  from  Controlling  Office)  ] 15.  SECURITY  CL  ASS.  (of  this  repQT%_, 


UNCLASSIFIED 


i7  Fl 


15a  DECLASSIFICATION  DOWNGRADING 
SCHEDULE 


16  distribution  STATEMENT  (of  th/a  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.  Distribution  statement  rof  the  sbsrract  entered  }n  Block  20.  if  different  from  Report) 


D D C 

HP' 

SEP  2 1977 

feLbiru  ni 


I 18  supplementary  notes 


Contains  information  presented  at  the  Symoosium  on  Air  Force  Applications  of 
Modern  Control  Theory,  July  1974  and  at  the  Modeling  and  Simulation  Confer- 
ence, April  1975. 


19.  KEY  WORDS  f Continue  on  reverse  aide  if  necaaaarv  ami  Identify  hy  block  number) 

Numerical  Integration 

Numerical  Solution  of  Differential  Equations 
Tunable  Integration 
Control  Theory 

Error  Analysis  

20  ABST  lu#  on  reverse  sidi  if  neceeeery'  end  ide  'ttfy  by  block  numbe") 


F'^equency  Response 
Root  Locus 
Bode  Plots 
Stabi 1 ity 


Tunable  integration  is  an  aoproach  to  the  numerical  solution  of  ordinary 
differential  equations  that  is  conceived,  developed,  and  employed  based  upon 
the  principles  of  linear  control  theory.  -<his  present  work  complements  the 
original  publications  on  tunable  integratiort'-by^mith,  and  it  also  reiter- 
ates and  extends  the  earlier  work  hy  this  author. The  properties  of  the 
zero-order-hold  tunable  integT'ator  and  the  control -theory  methods  used  to 
analyze  those  properties  are  equally  important  in  the  context  of  this  report. 


DD  . 1473  EOI'FION  OF  1 NOV  *5  'S  OBSOLF.TF 


UNCI.ASSIFIFD  ^ I J . 

9ECUBI'’Y  CL  ASS' F|C  F ’’ION  OF  -HIS  PAGE  rl*'h»r  Data  F.nttrad 


r 

UNCLASSIFIED 

1 

security  classification  of  ^ ,i\s  PAGE(»Vhefi  Date  Entered} 

20.  ABSTRACT 


Beginning  with  a discussion  of  numerical  error  from  a control -theory  perspec- 
tive, present  the  concept/of  an  Ideal  integrator  having  no  error  but  that 
due  to  finite  computer  word  length.  The  lack  of  adaptability  of  the  ideal 
integrator  motivates  the  tunable  integrator,  which  possesses  the  requisite 
flexibility.  The  majority  of  the  report  covers  the  analytic  development  of 
tunable  integration,  the  frequency  response  of  the  zero-order-hold  tunable 
integrator,  and  a root-locus  analysis  of  that  integrator  as  employed  in  a 
first-order,  linear  system.  Throughout  the  presentation,  emphasis  is  placed 
on  the  use  of  these  control -theory  techniques  to  determine  stability  and  to 
tune  the  integrator  for  improved  performance,  ^he  final  section  of  the  report 
stands  somewhat  apart  from  the  rest,  but  is  equaTTy  important  in  that  it 
discusses  the  effects  of  computer-code  implementation  upon  the  frequency 
response  of  a software  package.  The  results  of  the  studies  reported  on  here- 
in are  only  a start.  The  potential  promised  by  tunable  integration  deserves 
and  requires  further  investigation. 


unclassified 

security  Cl  ASSI  EIC  AT!0N  this  PACE'i^'hen  Enferetf) 


FOREWORD 


This  report  was  prepared  by  the  Applied  Mathematics  Division,  Directorate 
of  Aerospace-Mechanics  Sciences,  Frank  J.  Seiler  Research  Laboratory,  United 
States  Air  Force  Academy,  Colorado.  The  work  was  initiated  under  Project  Work 
Unit  No.  2304-F1-55,  "Tunable  Integration."  Major  Marc  L.  Sabin  was  the 
Project  Engineer-in-Charge. 

I would  like  to  give  particular  thanks  to  Mr.  J.  M.  Smith,  whose  comments 
and  counsel  were  invaluable  during  the  course  of  this  effort.  I would  also 
like  to  give  thanks  to  Capt  R.  P.  Fuchs  of  the  Department  of  Astronautics  and 
Computer  Science,  USAF  Academy,  for  his  many  suggestions  during  the  research 
and  for  his  critical  review  of  the  manuscript.  My  sincere  thanks  go  to 
Mrs.  D.  Weiss  for  her  efforts  in  typing  the  manuscript  and  preparing  it  for 
publication. 

This  technical  report  has  been  reviewed  and  is  approved  for  publ icat^on. 

''  ... 


...  - 

MARC  L.  SABIN,  Major,  U.SAF 
/'Deputy  Director 
Aerospace-Mechanics  Sciences 


BURTON  H.  HOLAOAY,  Lt  Col,  USAF  ' 
Di rector 

Aerospace-Mechanics  Sciences 


FOR  THE  COMMANDER  , 

TFfOMAS  M.  TOMASKOVIC 
Chief  Scientist 


Arc''""’! 

f'.'.r 

/ 

»J-:S 

ction 

nr'. 

£ .!  ' .r. 

'on  1 1 

!•'  ■■ 

n 

'' 

1 

t 

1 PY 

I i;y 1 1 0 1 

"k 

si'i.cr.'- 1 

^ i 

J 


TABLE  OF  CONTENTS 


Page 

Foreword  1 

List  of  Figures  3 

List  of  Tables  A 

List  of  Symbols  and  Acronyms  5 

Chapter  1:  Introduction  7 

1.1  The  Purpose  of  This  P.eport 7 

1.2  The  Importance  of  the  Subject 7 

1.3  Tunable  Integration  8 

1 .4  Overview  9 

Chapter  2:  A Question  of  Errors  10 

2.1  Some  Elementary  Preliminaries  10 

2.2  An  Impoi^tant  Aside 12 

2.3  Errors  From  the  Control  System  Point  of  View 13 

2.4  The  True  Integrator 14 

2.5  The  Ideal  Integrator 17 

Chapter  3:  Introduction  to  Tunable  Integration  22 

3.1  Basic  Formulation  22 

3.2  Development  of  the  Zero-Order-Hold  Tunable  Integrator  26 

3.3  An  Example  Aoplication  29 

3.4  The  Advantages  of  Tunable  Integration  35 

Chapter  4:  Frequency-Response  Characteristics  of  the  ZOH-TI  36 

4.1  Why  the  Frequency  Response  36 

4.2  Methods  of  Evaluation  36 

4.3  Frequency  Response  of  the  ZOH-TI  37 

4.4  Tuning  Based  Upon  the  Frequency  Response 41 

Chapter  5:  Root-Locus  and  Stability  Considerations  44 

5. 1 The  Approach  44 

5.2  Root  Locus  of  a First-Order  System  with  the  ZOH-TI  44 

5.3  Stability  Considerations  — 50 

Chapter  6:  The  Effect  of  Implementation  59 

6.1  Wny  Worry? 7 59 

6.2  Implementation  of  tx  + x = f 59 

6.3  Comparison  of  Frequency  Responses  64 

6.4  The  Conclusion  66 

Chapter  7:  Conclusions  and  Recommendations  68 

7.1  Summary  68 

7.2  Conclusions  69 

7.3  Recommendations  69 

Bibliography  71 


2 


LIST  OF  FIGURES 


Page 

1.  Quadratic  Fit  to  y II 

2.  Error  vs  Step  Size  12 

3.  The  Numerical  Integrator  as  a Control  System  14 

4.  True  Integrator  Pole  Location  15 

5.  True  Integrator  Frequency  Response  16 

6a.  Continuous  Compensation  22 

6b.  Discrete  Compensation  22 

7.  Stair-Step  Approximation  23 

8.  Zero-Order  Hold  Impulse  Response 24 

9.  ZOH-TI  26 

10.  ZOH-TI  Pole-Zero  Map - - 27 

11.  Second-Order  System  Block  Diagram  30 

12.  ZOH-TI  Solution  for  y = 0.0 31 

13.  ZOH-TI  Solution  for  y = 1-0 32 

14.  ZOH-TI  Solution  for  , = 0.5 33 

15.  ZOH-TI  Step-Response  Solution  34 

16.  ZOH-TI  Bode  Magnitude  Plot 40 

17  ^OH-II  Bode  Phase  Angle  Plot 41 

H-Tl  Step  Response  with  y = 0.52  43 

■irst-Order  System  Block  Diagram  44 

20.  ZOH-TI/Euler  Implementation  of  ix  + x = f 47 

21.  Rectangul ar/Euler  Implementation  of  tx  + x = f 48 

22.  Root-Locus  for  ZOH-TI/Euler  Implementation  50 

23.  Stability  Region  for  Tuning  52 

24.  True  Euler  Implementation  of  tx  + x = f 56 

25.  Mixed-Mode  Euler  Implementation  of  tx  + x = t 61 

26a.  True-Euler  Implementation  of  tx  + x = f 63 

26b.  Reduced  Block  Diagram  63 

27.  Comparison  of  FT-equency-Response  Characteristics  67 


3 


LIST  OF  TABLES 


Page 


1.  Ideal  Cosine  Integrator  Frequency  Response 18 

2.  Ideal  Cosine  Integrator  Performance 19 

3.  Ideal  Integrators  20 

4.  The  ZOH-TI  and  Classical  Integrators  29 

5.  True  Solution  vs  ZOH-TI/Euler 54 

6.  Frequency-Response  Comparison  56 


LIST  OF  SYMBOLS  AND  ACRONYMS 


Engl ish 


Arg  ■'*} 

C(s) 

E 

= 9f/3x 
G(ju)),  G*(j(ji)) 

G(s),  G(z) 

GH 

H 

K 

M 

Mjjb  = 20  log^Q  M 
p = ATy 
n = (y-I)/i 
Rc(s),  Ri(s) 

Rzoh^^) 

r = 

s 

T 

^CL 

TI 

t 

ii(t) 

!■*] 

ZOH-TI 
7 ' 


Phase  angle  of  the  complex  quantity  I*} 

Compensator  transfer  function 
Total  local  error 
Local  roundoff  error 
Local  truncation  error 
Partial  derivative 

Frequency-response  and  sampled  frequency-response  transfer 
functions,  respectively 

Laplace  and  sampled-data  transfer  functions,  respectively 
Open-loop  transfer  function 
Feedback  transfer  function 
Forward  path  gain 

Magnitude  of  frequency-response  transfer  function 

Magnitude  in  db 

ZOH-TI  gain 

ZOH-Tl  zero  location 

Compensation  and  integration  reconsrructor  transfer  functions, 
resoecti vely 

Zero-o'der  hold  reconstructor  transfer  function 

Sampling  ratio 

Laplace  transform  variable 

Integration  step  size 

Closed-loop  transfer  function 

Tunable  integration  or  tunable  integrator 

Time 

Unit  step  function 
Z-trans'^orm  of  (*} 

Zero-order-hold  tunable  integrator 
Z-trar,sform  variable 


5 


Greek 

] 

Y 

TI  phase-compensation  parameter 

i 

A 

Forward  difference  operator 

Damping  coefficient 

X 

TI  gain-compensation  parameter 

: 

System  time  constant 

4> 

Argument  of  frequency-response  transfer  function 

ik(z) 

Characteristic  polynomial 

U) 

Frequency 

■ 1 

‘^d 

Design  frequency  or  damped  natural  frequency 

Natural  frequency 

uig  = 2tt/T 

Sampling  frequency 

Subscripts 

£l 

1; 

c 

Computed  value 

k 

kth  increment  or  value  at  t = kT 

mm 

Mi xed-Mode-Eul er  computed  value 

• 

t 

True  value 

te 

True-Euler  computed  value 

Other 

€ 

Inclusion;  e.g. , a € [b,c] 

1*1 

Magnitude;  e.g.,  M = |G(jo))| 

. 1 

Hr 

Derivative;  e.g.,  y = dy/dt 

1 



6 

] 

■A 

( 

chapter  1 

Introduction 

"The  Purpose  or  Computing  is  Insight,  Not  Numbers"-^ 

1.1  The  Purpose  of  This  Report 

The  principal  objective  of  this  report  is  to  generate  some  of  the  insight 
that  Dr.  Hamming  succinctly  stated  to  be  the  purpose  of  all  computation. 

Specifically,  I will  present  a novel  technique  of  numerical  integration  for 
the  solution  of  ordinary,  first-order  differential  equations  occurring  in 
ini ti al -val ue  problems;  and  I will  analyze  that  method  from  a control - 
theoretic  viewpoint  that  -is  applicable  to  the  analysis  of  any  numerical 
algorithm. 

Intended  readers  include  control  engineers,  who  can  greatly  extend  the 
efforts  discussed  herein,  and  numerical  analysts,  who  can  use  these  concepts 
along  with  their  own  effective  methods  to  further  enhance  our  ability  to 
numerically  solve  differential  equations. 

A third  audience  to  whom  this  repo*"!  is  addressed  is  comprised  of  the 
faculty  and  cadets  of  the  United  States  Air  Force  Academy  who  are  involved 
in  teaching  and  studying  numerical  methods.  Specific  aopl icabi 1 ity  is 
intended  for  courses  in  the  Departments  of  Astronautics  and  Computer  Science 
and  of  Mathematics.  To  be  useful  to  these  people,  the  style  of  this  repo*"t 
is  tutorial,  and  much  is  said  that  will  be  obvious  and  elementary  to  the 
reader  experienced  in  the  fields  of  control  theory  and/or  numerical  analysis. 

1.2  The  Importance  of  the  Subject 

The  United  States  Air  Force  could  not  oerform  its  mission,  it  could  not 
"fly  and  fight,"  if  it  did  not  possess  the  means  to  numerically  compute  the 
solutions  to  differential  equations.  The  motions  of  a hand  grenade  thrown 
by  a soldier,  a bullet  ♦^ired  from  a rifle,  an  aircraft  on  a strafing  pass,  a 
smart-bomb  following  a laser  designator,  a ballistic  missile  on  an  inter- 
continental trajectory,  or  a satellite  orbiting  the  earth  are  all  governed  by 
differential  equations. 

Differential  equations  describe  more  than  the  motion  of  objects.  The 
Riccati  equation  for  optimal  contro’  is  a di "'ferential  equation,  equations  of 
chemical  kinetics  are  differential  equations,  the  covariance  matrix  of  a 
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Kalman  filter  propagates  accoraing  to  a differential  equation,  and  electric 
circuits  are  described  by  differential  equations;  the  list  is  endless. 

In  most  realistic  situations,  analytic  solutions  cannot  be  found.  The 

need  for  numerical  solutions  is  self-evident;  they  are  needed  to  accomplish 

operational  missions,  they  are  needed  for  simulation  (especially  important 

v/ith  the  increasing  amount  of  training  via  simulation),  and  they  are  needed 

for  design  and  analysis.  Different  situations,  however,  require  different 

approaches  to  their  solutions.  The  formalized  differential-equation-solver 

packages  such  as  DVDO,  GEAR,  etc.,  are  excellent  for  highly  accurate 

solutions  on  large,  ground-based  systems.  For  the  small  airborne  computer, 

simpler,  less  generalized  approaches  are  required,  such  as  the  technique 

( 3 1 

employed  in  the  ballistic  trajectory  algorithm  of  Duke,  et.  al.'^  ' 

1.3  Tunable  Integration 

A number  of  years  ago,  J.  M.  Smith  thought  it  reasonable  that  since  the 

digital  computer  is  a discrete  operating  device,  the  techniques  of  sampled- 

data  system  analysis  and  digital  filtering  could  effectively  be  used  to 

design  methods  for  numerical  integration  on  the  digital  computer.  Thus  was 

born  the  concept  of  tunable  integration  as  it  is  presented  herein.  Smith 

141 

has  used  the  technique  with  considerable  success.^  ' Only  little,  however, 
has  been  formally  written  concerning  the  method^ , and  this 
report  is  intended  to  increase  its  familiarity  within  the  Air  Force  com- 
munity. 

The  technique  is  based  upon  the  synthesis  of  a discrete  approximation 
to  continuous  integration.  The  result  is  a difference  equation  in  which 
adjustable  parameters  orovide  the  capability  to  directly  control  the  charac- 
teristics (pole  location,  frequency  response)  of  the  integrator.  This 
capability  makes  tunable  integration  a technique  that  proffers  tremendous 
benefits  in  terms  of  ease  and  flexibility  of  use,  large  regions  of  stability, 
controllable  accuracy  (via  means  other  than  variable  order  or  step  size), 
and  simplicity.  High-order  integrators  of  this  type  have  been  developed  by 
Smith^^^ ’ , but  this  author  feels  that  more  significant  benefits  accrue  to 
the  Air  Force  from  the  use  of  the  low-order  form  in  appropriate  problems  such 
as  airborne  fire  control.  Thus,  th’s  report  will  restrict  its  attention  to 
the  low-order  form  of  tunable  integration. 
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The  organization  of  this  report  is  based  on  using  tunable  integrafion  as 
a vehicle  for  demonstrating  the  application  of  control  theory  in  analyzing 
methods  of  numerical  integration.  The  integ>^ation  algorithm  itself  and  the 
method  of  analysis  are  equally  important  parts  of  this  reoort. 

1.4.1  A ef  Outl  i ne 

Chapter  2 presents  a general  discussion  of  errors  and  how  to  define 
error  criteria  in  terms  of  the  pole  location  and  frequency  response.  The 
concept  of  "ideal"  integrators  is  used  to  demonstrate  these  error  criteria 
and  to  motivate  the  tunable  integrator.  Though  the  development  also  appears 
elsewhere ‘ ' , the  complete  formulation  of  the  zero-order-hold  tun.-ble 

integrator  is  presented  in  detail  in  Chapter  3.  An  example  problem  i^  solved 
using  that  integrator  and  clearly  demonstrates  the  tuning  property.  Chapter 
4 presents  the  frequency-response  characteristics  of  the  integrator.  In 
Chapter  5 a root-locus  analysis  is  performed  and  stability  is  evaluateo. 

Also  in  that  chapter  is  a short  demonstration  of  how  the  frequency- response 
and  root-locus  analyses  can  be  used  to  initially  tune  an  integrator.  The 
importance  of  program  implementation  is  discussed  in  Chapter  6,  wherein  the 
difference  in  properties  between  two  program  codes  using  the  same  basic 
method,  but  programmed  differently,  is  shown.  Shampine,  et.  al.,'  ‘ have 
noted  the  significance  of  differing  implementation  upon  the  results  obtained 
with  GEAR  and  DIFSUB,  two  packages  using  the  same  fundamental  methods. 

1.4.2  Background  Knowledge 

In  o'^der  to  understand  this  report,  the  reader  should  be  familiar  with 
sorif  basic  concepts  of  linear  and  sampled-data  control  systems,  and  of 
numerical  analysis.  Specifically,  the  reader  should  be  somewhat  familiar 
with  tue  following;  Laplace  transforms,  block  diagrams,  transfer  functions, 
root-locuS  analysis,  frequency-response  analysis,  control -system  stability, 
phase  and  gain  compensation,  sampling  of  continuous  signals,  the  sampling 
theorem,  7-transformations , difference  euuations,  polynomial  interpolation, 
and  classical  methods  of  numerical  integration.  The  true  integrator  is 
continuous  and  the  numerical  integrator  is  discrete.  Thus,  familiarity  with 
concepts  in  both  continuous  and  discrete  control  theory  is  important. 


CHAPTER  2 

A Question  of  Errors 


In  his  "Essay  on  Numerical  Methods" Hamming  cites  five  main  ideas 
germane  to  the  study  of  numerical  methods.  The  first  of  these  is  quoted  at 
the  opening  of  Chapter  1.  Two  other  ideas,  or  concepts,  are  truncation  and 
roundoff  error,  the  fundamental  error  sources  of  any  numerical  method  to  be 
implemented  on  a digital  machine.  Perhaps  the  most  important  consideration 
for  any  numerical  metnod  is  its  error  characteristics.  The  methods  of  control 
theory  provide  a different  perspective  on  the  errors  inherent  in  a numerical 
algorithm.  Tunable  integration  is  a method  of  numerical  integration  that  is 
designed  to  take  advantage  of  that  perspective  and  to  thereby  provide  excel- 
lent error  characteristics.  In  this  chapter  we  shall  look  at  the  fundamental 
question  of  errors. 

2.1  Some  Elementary  Preliminaries 

Roundoff  error  arises  from  the  finite  length  of  the  "word"  in  a digital 
computer.  Hence,  any  number  that  cannot  be  represented  by  a finite  number  of 
binary  bits  (most  numbers)  must  be  rounded  off  by  the  computer.  The  problem 
is  unavoidable,  and,  as  Hamming  suggests,  one  can  only  do  one's  best  to  mini- 
mize its  effects.  A principal  problem  with  'roundoff  is  that  it  tends  to 
worsen  as  the  number  of  computations  increase  or  as  computer  word  length 
decreases.  Thus,  roundoff  can  be  especially  severe  when  using  a mini  or 
micro  computer  to  solve  a differential  equation  over  a relatively  large 
(compared  to  the  integration  step  size)  interval. 

Truncation  error,  on  the  other  hand,  tends  to  require  smaller  integra- 
tion steps,  since  it  arises  from  the  necessity  to  use  finite  polynomial 
expansions  to  represent  arbitrary  functions.  In  a typical  classical  approach 
to  the  development  of  numerical  integrators,  the  integrand  is  expanded 
polynomially  and  a truncated  series  from  that  polynomial  is  analytically 
integrated  over  a specified  interval  to  provide  the  desired  integration 
formula. 
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As  an  example,  let  us  derive  Simpson's  integration  formula.’  ' Newton's 
forward  interpolation  formula  can  be  written  as 


■W 


10 


Pk  ' £ (?) 


n - k 


= yo  kAy^  + ‘5  k(k-l)A2y^  + 


(2-1) 


where  Pj^  is  the  polynomial  expansion  of  y to  n forward  differences, 

Ay^  = y^^^  - y^.  Let  n = 2 so  that  we  are  approximating  y by  a quadratic  in 
k over  the  three  point  interval  y^,  y^,  y^  as  shown  in  Figure  1 below. 


Figure  1.  Quadratic  Fit  to  y 


The  integration  over  the  interval  gives  us 


f ';(t)dt  = y^-y^  = if  P^c 


(2-2) 


where  we  note  that  dt  = Tdk.  Substituting  for  P|^  from  Eq  (2-1)  and  then 
expanding  the  forward  differences,  we  obtain 


/•2  _ . 

y?  " yo^'j  k(k-l)A2y^]  dk  = y^  a - (y^  + Ay^+yJ 


(2-3) 


Equation  (2-3)  is  Simpson's  rule,  which  has  the  general  form 

V2  ■ I * ’Vl  * >'n> 

To  obtain  it,  we  approximated  the  integrand  y by  a polynomial  truncated  to 
the  second  order  (in  k)  and  then  analytically  integrated  to  obtain  the 
desired  result. 

Truncation  error  arises  from  all  of  the  neglected  terms  in  the  expansion 
of  P|^.  It  represents  the  error  in  approximating  the  true  integrand  by  the 
selected  polynomial  over  the  interval  of  integration.  Clearly,  the  larger 
the  interval  over  which  a given  order  polynomial  is  made  to  fit  an  arbitrary 
function,  the  larger  the  errors  will  tend  to  be,  and  the  more  significant 
will  be  the  truncated  terms. 

The  contradictory  requirements  of  roundoff  and  truncation  error  are 
shown  in  Figure  2. 


ERROR 


NUMBER  OF  COMPUTATIONS  FOR  FIXED  INTERVAL 
Figure  2.  Error  vs  Step  Size 


2.2  An  Important  Aside 

We  have  just  derived  the  formula  for  Simpson's  Rule  by  truncating  the 
Newton  interpolation  formula.  Unless  the  function  y itself  goes  to  infinity 
over  the  interval  [t^,  t^,],  the  error  in  the  approximation  will  be  finite. 
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and  the  error  over  a finite  nuirber  of  intervals  will  also  be  finite.  The 
integration  of  y,  however,  is  here  a quadrature;  a given  y,  explicitly  defined 
as  a function  of  t,  has  been  approximated  over  a known  interval.  In  the  more 
general  problem  with  dynamics,  y will  depend  on  y,  and  errors  in  the  computa- 
tion of  y will  be  fed  back  ^o  propagate  and  perhaps  grow  without  bound.  In 
such  a situation,  an  unbounded  error  can  arise  in  a finite  interval. 

The  question  of  such  growth  in  errors  is  a question  of  stability,  and  we 
will  see  later  how  the  control  system  perspective  provides  a very  convenient 
way  to  examine  the  stability  of  a numerical  integrator.  Note  for  now  that 
many  integrators,  of  which  Simpson's  is  one,  may  be  acceptable  for  quad'^ature 
problems  but  are  patently  unstable  for  the  dynamic  problem. 


2.3  Errors  From  the  Control  System  Point  of  View 

From  the  discussion  above,  we  can  see  how  the  classical  approach  leads 
us  to  the  use  of  the  integration  step  size  or  the  order  of  the  integrator 
(i.e.,  the  degree  of  truncation)  to  the  cc  .trol  of  errors.  Means  are  avail- 
able to  evaluate  bounded  magritudes  for  the  errors,  and  can  be  found  in  many 
references  f 1 J , [ 10] , 1 1 1 | . When  one  uses  variable-step-size  or  variable-order 
algorithms  to  control  those  bounds,  one  pays  the  penalty  of  problems  in 
increased  complexity  and  increased  computational  burden  as  the  order  increases. 
Tunable  integration  and  the  control -theory  analysis  provide  an  alternative 
approach  to  error  control,  one  in  which  neither  step  size  nor  order  are 
varied.  The  parameters  that  tune  the  integrator  control  the  error.  Avoided 
are  the  restart  problems  of  changing  step  size  (for  the  multistep  methods) 
and  the  complexity  of  higher  order. 

To  discuss  errors  from  the  control -system  perspective,  we  must  look  at 
the  integrator  itself  as  a control  system.  If  we  define  a control  system  as 
a collection  of  components  and/or  algorithms  designed  to  take  a given,  known 
or  uncertain,  input  signal  and  from  it  produce  an  output  with  certain  desired 
properties,  then  the  function  of  an  integrator  is  clear:  the  numerical  inte- 
grator is  a 1 inear-transfer- function  control  system  that  takes  an  input 
integrand  and  produces  as  its  output  the  integral  of  the  input.  Schematically 
we  can  portray  this  function  as  in  Figure  3. 
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Figure  3,  The  Numerical  Integrator  as  a Control  System 

The  true  integrator  would  produce  the  true  integral  at  the  output.  How 
far  from  the  true  integrator  the  numerical  integrator  is  can  be  analyzed  in 
terms  of  the  transfer  function  (using  the  Laplace  variable  s): 

G(s)  = 

There  are  two  fundamental  properties  of  G(s),  the  analyses  of  which  comprise 
the  heart  of  this  report.  These  properties  are  1)  the  locations  of  its  poles 
and  zeros,  and  2)  the  magnitude  and  phase  characteristics  of  its  frequency 
response.  How  these  properties  differ  from  these  of  the  true  integrator  is  a 
direct  measure  of  the  performance  of  our  numerical  integrator. 

2.4  The  True  Integrator 

In  the  Laplace  domain  the  transfer  function  of  the  true  integrator  is 

Gj(s)  = i (2-5) 

This  single  pole  at  the  origin  maps  into  a single  pole  at  (1,0)  in  the 
z-plane,  where  z = e^^.  Thus,  pole-zero  maps  of  the  true  integrator  in  the 
s-  and  z-planes  are  simply  as  shown  in  Figure  4. 
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Fiaure  4. 


True  Integrator  Pole  Location' 


We  should  like  any  numerical  algorithm  to  have  a pole  close  to  that  of  the 
true  integrator.  We  should  also  like  as  few  additional  poles  in  the  algorithm 
as  possible,  for  additional  poles  result  in  additional  lag  which  tends  to  be 
destabilizing,  and  they  add  to  the  complexity  of  the  algorithm. 

A second  basic  property  of  the  true  integrator  is  its  frequency  response 
magnitude  and  phase.  Substituting  ju  for  s ■’’n  Eq  (2-5),  we  have 


G^(jal) 


The  magnitude  and  phase  are 


(2-6) 


It  is  interesting  to  note  that  simply  taking  the  z-transform  of  ^ would 
give  G^(z)  = Z j-j  = The  corresponding  difference  equation  is  = 

+ ^n+1 ’ which  implies  a unit  integration  step  in  order  to  have  any 
practical  meaning. 


15 


M = |Gjjco)|  = I 

= 20  log,„M  = -20  log  uj 
db  ^10  ^10 

<>  = Arg  {G^(ju)}  = -90° 

The  corresponding  Bode  diagrams  are  as  shown  in  Figure  fj. 


(2-7) 

(2-8) 


Figure  5.  True  Integrator  Frequency  Response 

Any  numerical  algorithm  should  have  magnitude  and  phase  characteristics 
closely  approximating  the  20  db/decade  gain  slope  and  90°  of  phase  lag,  at 
least  over  the  frequency  range  of  interest. 

These  two  fundamental  properties  of  pole  location  and  frequency  response 
provide  a control  system  perspective  for  the  analysis  of  numerical  integration 
techniques  that  will  be  used  throughout  the  remainder  of  this  report.  Varia-  | 

tions  in  pole  location,  additional  poles  and  zeros,  variation  of  M from  1/u, 
and  variation  of  from  -90°  are  all  direct  measures  of  the  error  in  the 
numerical  integration  algorithm.  We  will  later  see  a clear  tie  between 


The  use  of  these  criteria 
simulation. 


is  not  unfamiliar  to  the  engineer  involved  in 
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phase  and  gain  distortion  and  the  classical  concept  of  truncation  error. 

This  will  arise  from  the  use  of  holding  circuits  to  reconstruct  sampled 
signals.  The  reconstructed  signals  are  simply  polynomial  approximations  to 
their  continuous  predecessors. 

2.5  The  "Ideal"  Integrato'"’^ 

Having  defined  what  will  be  our  control  system  perspective,  we  can 
derive  a set  of  ideal  integrators  that  have  the  exact  properties  desired,  but 
only  for  the  design  input.  These  ideal  integrators  take  specified  inputs  and 

4- 

transform  them  into  the  corresponding  known  integrals  for  their  output. 

For  example,  consideran  ideal  cosine  integrator.  The  z-transforms  of  input 
and  output  are 


z (z  - cos  (iJ^T) 

Z {cos  W jt  / “ T \ T 

d z'  - (2  cos  (jj^T)  z+1 


Z Csin  oj^t)  = 


(sin  wjT)  z 


z-  - (2  cos  wjT)  z+1 


where  is  the  design  frequency.  The  ideal  transfer  function  is 


Gi(z) 


sin  o',^t 




Wj  Z{C0S  oj^tl  ' Xz  - cos  WjjT) 


(2-9) 


The  Bode  magnitude  and  phase  are  found  by  substituting  for  z in  Eg  (2-9) 

G'-(ju>)  = - 


sin  tj^T 


This  section  is  based  on  ideas  proposed  and  initially  developed  by  Ronald  t. 
Janosko. 

' Whis  concept  is  not  too  dissimilar  from  the  ideas  presented  by  Fowler. 

Note  that  we  have  the  sampled  frequency  response  G*(jj)  and  not  the 
continuous-system  response  G(ju). 


(2-10) 


sin 


[1  + cos^  - 2 cos  aiT  cos  u^T] 


. -k 


4)  = -tan 


-1  ) sin  (oT  ( 
I cos  u)T  - cos  (D  ,T  \ 


At  the  design  frequency  these  relations  reduce  to 


Thus,  at  the  design  frequency  we  have  the  desired  frequency-response  charac- 
teristics, independent  of  the  time  step  selected.  There  is  also  only  one 
pole,  though  its  location  will  depend  upon  the  value  of  u^T. 

What  happens,  however,  if  we  use  this  ideal  integrator  at  other  than  the 
design  frequency?  If  we  let  = ir/6  (the  sampling  ratio  is  thus  1/12;  i.e., 
“d^“s  ” “(j/(2’'/T)  = 1/12),  Eqs  (2-9)  yield  the  data  shown  in  Table  1. 

Table  1:  Ideal  Cosine  Integrator  Frequency  Response 


u.d-M 

U)  • M 

4>(deg) 

0.5 

1.80 

0.90 

- 69 

0.7 

1.37 

0.96 

- 79 

1.0 

1 .0 

1 .0 

- 90 

1.3 

0.79 

1 .02 

- 98 

1.5 

0.69 

1.04 

-103 

Both  the  magnitude  and  phase  characteristics  degrade  as  we  move  away  from 
Note  that  both  products,  . M and  w • M,  drift  away  from  unity,  which  the 
true  integrator  produces  at  all  frequencies. 

We  noted  above  that  the  frequency  response  at  the  design  frequency  was 
independent  of  T.  Concomitantly,  the  accuracy  is  independent  of  T;  the 
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ideal  integrator,  not  surprisingly,  gives  exact  results  regardless  of  step 
size.  This  can  be  seen  by  writing  the  difference  equation  for  the  integrator. 
From  Eq  (2-9)  that  difference  equation  is 


where 


X = ax  + bx 
n+1  n n 


a = cos  cj  ,T 
d 


b = — sin  0)  ,T 
“d 


Given  the  initial  conditions  x = cos  nuj.T  and  x = — sin  nw  ,T,  we  algebraically 

1 ^ 

obtain  the  result  x = — sin  [(n+l)ai.T],  independent  of  T.  Using  the  same 

example  as  was  used  to  generate  Table  1,  we  can  perform  the  numerical  integra- 
tion and  see  the  actual  errors.  Let  = 1 and  T = /6  so  that  oi^T  = n/6  as 

before.  The  data  in  Table  2 show  ten  integration  steps  each  for  w = 0.5, 
to  = 1.0  = to^ , and  ui  - 1.5. 


Table  2:  Ideal  Cosine  Integrator  Per-*"o>'mance 


t 

1.5 

1.0 

1.5 

■v' 

X * 
(. 

V 

sr  1 

"t 

X 

c 

n 

n 

c 

0 

0 

n 

n.'.n 

0.513 

0.5'iO 

0 . 500 

0.555 

0.5(71 

0.500 

1.  17 

1.000 

).916 

0.056 

O.l.fG 

0.66  7 

O./:' 

1.S71 

i.'ir. 

1 •O.-'O 

1 000 

1.  :'J0 

0.771 

( .•  M 

:.i  ya 

1.7J? 

l.-ilO 

0.'66 

0.:'55 

0.000 

0.  ’6 

c'.nin 

1 .93? 

1 .7/5 

0.5':) 

0.509 

-0.771 

-U..  ^5 

J.U.2 

?.000 

1 3 

0.330 

O.Oj" 

-.3.667 

J.f'lio 

1.93? 

1.219 

-.'1.5C9 

-0.'  '0 

-0.771 

-0.  . , 

■1.  ro 

1.732 

-0.765 

-C.Otf 

0.000 

-0.1  11 

7.71? 

1 . -1 1 .j 

0.' .5; 

-1.00J 

-1 .000 

0.471 

4.410 

5.27(> 

1 . 0'lO 

0.  1.'5 

-0.''5/ 

-0.366 

0.667 

•1.  /;(  1 

Notf.:  1.  ‘ '-sin  -1.^;  IruP  intpyiotcr  result 

^ c 'Sine  intcyrator  result 

sir  .jT 

X (t  ) -os  ..  t ; a ‘ cos  ...i  : b = 

c !i  n Vi  .1 
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The  relation  between  the  data  in  the  two  tables  is  clearly  seen  in  the  case 
of  tjj  = 0.5.  The  magnitude  is  too  small  and  the  phase  lag  is  too  little  in  the 
first  table.  The  peak  seen  in  Table  2 is  less  than  the  true  peak  and  occurs 
before  it  should,  just  as  the  frequency  response  data  indicated.  For  w = 1.5 
the  1 arger-than-ideal  magnitude  is  evident,  but  the  excessive  phase  lag  cannot 
be  seen  due  to  the  spacing  of  the  data  points. 

Some  additional  ideal  integrators  are  shown  in  Table  3. 

Table  3:  Ideal  Integrators 


Input 

G(z) 

u(t) 

T/(z-l) 

t 

e^t 

1/a  (e^^-1) 
z-1 

sin  u)^t 

1 - cos  o^T 

0)^  sin  w^T  ‘ z-1 

cos  m^t 

s i n (i!^T  1 

01^  z - cos  w^T 

*Tustin's  Transformation 


From  this  elementary  analysis,  admittedly  using  an  example  of  quadrature 
to  make  our  points,  it  is  clear  how  one  can  apply  the  methods  of  the  control 
engineer  to  the  analysis  of  a numerical  integrator.  We  have  seen  how  the 
ideal  integrator  may  be  derived  from  the  perspective  o-f"  a linear  transfer 
function.  Its  range  of  usefulness  is,  however,  very  limited.  If  we  define  a 
sampling  ratio  r as  follows 

r = Wg  = 2t/T  (2-11) 


. 

I 

i 


20 


then,  for  the  example  given  earlier,  w = 0.5  corresponded  to  r = 0.04,  ou  = 1 
to  r = 0.08,  and  w = 1.5  to  r = 0.13.  We  shall  see  ■'n  Chapter  3 that  this  is 
indeed  a very  limited  range  of  the  sampling  ratio.  If  one  were  to  design  an 
integrator  from  the  ground  up  using  a control -theory  perspective,  a prime 
consideration  would  be  making  the  useful  range  of  r as  large  as  possible. 

This  is  precisely  what  has  been  done  in  the  development  of  the  tunable  inte- 
grator which  we  will  discuss  in  depth  beginning  with  the  next  chapter. 


21 


CHAPTER  3 

Introduction  to  Tunable  Integration 

This  chapter  begins  the  heart  of  the  report:  the  analysis  of  ti  lable 

integration  (TI).  In  it  we  will  briefly  review  the  formulation  of  TI,  which 
is  the  brainchild  of  Jon  M.  Smi th. We  will  derive  the  zero-order-hold 
TI  (ZOH-TI)  formula  and  will  demonstrate  its  application  to  a second-order 
differential  equation.  The  following  two  chapters  will  present  frequency- 
response  and  root-locus  analyses,  respectively. 

3.1  Basic  Formul atiori 

In  Smith's  words,  the  approach  to  developing  the  TI  is  one  of  "Synthesizing 
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a discrete  approximation  to  continuous  integration.  . ' To  do  this,  four 

basic  components  are  required;  1)  samplers  to  represent  the  discretization  of 
the  entire  process  of  digital  computation;  2)  a reconstructor  to  provide  a 
continuous  signal  from  the  sampled  input;  3)  a compensator  to  enable  control 
of  the  distortion  introduced  by  both  the  samoling  and  reconstruction  processes; 
and  4)  an  integrator  to  integrate  the  compensated  signal.  Given  these  four 
components,  there  are  two  ways  to  construct  the  discrete  integrator.  The 
first  is  to  use  continuous  compensation  as  shown  in  Figure  6a,  and  the  second 
is  to  use  discrete  compensation  as  shown  in  Figure  6b. 


Figure  6a.  Continuous  Compensation 


Figure  6b.  Discrete  Compensation 
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Two  families  of  integrators  have  been  generated  by  Smith  depending  upon 
the  type  of  compensation.  In  general,  -For  the  same  compensator,  the  continuous- 
compensation  integrators  are  less  complex.  Formulas  for  both  types  are 
presented  in  [5J.  Reference  [6]  develops  a set  of  continuous-compensation 
foniiulas  that  are  less  complicated  than  Smith's.  In  this  report  I will 
discuss  almost  exclusively  the  ZOH-TI.  Let  us  now  look  in  detail  at  the  recon- 
structor and  compensator. 

3.1.1  The  Reconstructor 

As  stated,  the  reconstructor  provides  a continuous  signal  from  the  samples 
it  is  provided.  The  complexity  of  the  reconstructed  signal,  and  its  amplitude 
and  phase  distor'tion,  will  be  a function  of  the  order  of  the  reconstructor. 

That  order  is  essentially  the  order  of  the  polynomial  that  is  fit  to  the  sample 
points.  Since  we  are  here  fitting  polynomials  to  data  points,  this  process  is 
related  to  the  process  typically  followed  in  classical  developments  of  numeri- 
cal integrators.  The  difference  arises  from  our  control -system  perspective 
and  the  use  of  a compensator  which  allows  us  to  control  the  errors  induced  by 
the  polynomial  approximation. 

The  zero-order-hold  reconstructor  provides  a stair-step  approximation  to 
a continuous  signal  as  shown  in  Figure  7. 


Figure  7.  Stair-Stec  Aooroximation 
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If  we  examine  the  response  to  a single  sample,  we  see  that  it  is  simply  a 
rectangular  pulse  of  magnitude  equal  to  the  sample  value  and  of  duration  T. 
Since  the  sample  is  ideally  an  impulse,  that  rectangle  is  the  impulse  response 
of  the  zero-order  hold  as  seen  in  Figure  8 (for  a unit  impulse). 


'’ZOH 


The  figure  also  shows  that  the  impulse  response  ^ posi- 

tive step  function  at  !;  = 0 and  a negative  step  at  t - T : 

rzoH(t)  = u(t)  - u(t-T)  (3-1 ) 

Taking  the  Laplace  transform  of  Eq  (3-1)  gives  the  standard  expression  for  a 
zero-order-hold  transfe'^  function: 


^ZOH 


(3-2) 


We  could  take  the  z-transform  of  R and  find  that  helps 

little  since  the  transformation  that  must  be  taken  for  continuous-compensation 
TI  is  ZfR^C  1}. 

Other  reconstructors  that  have  been  used  to  develop  Tl  formulas  are  first- 
and  second-order  and  triangular  holds,  the  last  not  being  physically  realize- 
able  but  still  being  mathematical ly  usei^ul.  While  capable  of  improving 
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accuracy,  high-order  holds  add  complexity  and  additional  poles,  both  being 
cont'^ary  to  the  objectives  of  simplicity  and  minimization  of  extraneous  poles. 
This  is  one  of  the  primary  reasons  that  this  report  deals  primarily  with  the 
ZOH-TI. 

3.1.2  The  Compensator 

The  function  of  the  compensator  is  to  enable  control  of  the  distortion 
introduced  by  the  sampling  and  reconstruction  processes.  Rather  than  utilize 
the  standard  form  of  lead-lag  compensators.  Smith  employs  the  complex  exponen- 
tial form 

C(s)  = xe"^^^  (3-3) 

Tfe  'wn  parameters  X arid  7 will  be  seen  later  to  be  the  gain  and  phase  compen- 
sation parameters  of  the  TI.  If  we  let  s = a + jj,  then 

]C(j.o)l  = Ae^'’^ 

Arg  xc(ju)}  = v^T 

These  nelations  show  us  that  while  > pnimarily  affects  phase  compensation,  it 
also  affects  the  gain.  We  will  see  this  explicitly  in  the  difference  equation 
fon  the  ZOH-T I . 

In  employing  the  compensator.  Smith  absorbs  e'^^  ’nto  the  product  R^C  j, 
transforms  back  to  the  lime  domain  and  then  transforms  the  result  to  the  Z- 
domain.  I have,  instead,  expanded  the  exponential  and  truncated  after  the 
second  term  (admittedly  creating  more  truncation  error,  but  the  control 
perspective  allows  me  to  monitor  and  counter  the  deliterious  effects).  It  is 
in  'lis  r.arner  that  the  continuous-compensation  integrators  of  [6)  were 
developed.  The  compensator  used  for  the  remainder  of  this  report  is  thus 

C(s)  - X(1  + >sT)  (3-4) 

It  is  interesting  that  noth  Smith's  approach  and  my  approach  give  the  same 
result  for  the  70H-TI.  The  discrete  compensator  used  by  Smith  is  developed 
using  a triangular  hold  and  the  compensator  of  Eq  (3-4).  When  used  as  in  [5| 


with  (see  Figure  6b)  being  a zero-order  hold,  the  result  is  again  the  same 
form  as  the  continuously  compensated  ZOH-TI. 

3.2  Development  of  the  Zero-Order-Hold  Tunable  Integrator 

3.2.1  The  Transfer  Function 

Using  continuous  compensation,  a zero-order-hold  reconstructor  and  the 
truncated  polynomial  compensator,  we  can  redraw  Figure  6a  as  follows: 


Figure  9,  ZOH-TI 


Note  that  the  input  X^(s)  is  sampled  to  provide  X,*(s),  and  the  continuous 
output  X,(s)  is  sampled  to  provide  X.*(s).  The  transfer  function  of  the 

c. 

continuous  path  between  the  samplers  is  simply  the  product  of  the  three 
elements  in  the  path; 


X(1  +yST) 


1 

s 


By  taking  the  z transformation  of  this,  we  obtain  the  total  transfer  function 
between  the  input  and  output  of  this  discrete  integrator.  Initially  we  can 
wri  te 


GU) 


z-1  7 < 1 + ysT  ) 
^ I s‘  1 


From  a set  of  transform  tables  such  as  those  in  [10]  and  [11],  we  have  the 
resul t 


''Though  no  sampler  exists  between  the  reconstructor  and  compensator,  the 

mathematics  of  z-transforms  allows  us  to  factor  out  the  term  (1-e  the 

trans^n relation  '•.'•lich  i<;  (z-l)/z. 


G(z)  = 


> 


(3-5) 


= l(z^ 

z-1 


where 


(3-6) 


Equation  (3-5)  is  the  transfer  function  of  the  ZOH-TI.  Equations  (3-6) 
reiterate  the  earlier  point  that  >•  affects  both  the  phase  and  gain  of  the 
integrator.  More  will  be  said  about  the  frequency  response  and  root-locus 
characteristics  of  this  integrator  in  the  next  two  chapters,  but  for  the 
moment  we  can  note  that  the  pole-zero  map  of  the  basic  integrator  (an  open- 
loop  device)  is  as  shown  in  Figure  10. 


Im{zl 


Re(z) 


Figure  10.  ZOH-TT  Pole-Zero  Map 


There  is  a pole  at  (1,0)  matching  the  true  integrator  shown  in  Figure  4. 

There  is  also  a zero  whose  location  depends  on  > according  to  the  relation 

for  q in  Eq  (3-6).  For  y 2.  ^^6  zero  is  always  within  or  on  the  unit  circle, 

while  for  y < 'i  the  zero  will  always  be  outside  the  unit  circle  (to  the  left 
for  Y > 0 or  to  the  right  for  y < 0,  though  the  desirability  of  ever  using 

Y ^ 0 is  very  doubtful  due  to  the  excessive  lag  that  would  be  induced).  If  a 

loop  were  placed  around  the  integrator  to  form  a simple,  first-order  system, 
then  increasing  the  forward  path  gain,  p,  by  increasing  T would  drive  the  pole 
to  the  zero.  For  y the  pole  would  never  exit  the  unit  circle  and  the 
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system  would  never  go  unstable  (it  would  be  inaccurate,  however),  while  for 
y > h’>  sufficiently  large  integration  steps  would  result  in  instability.  This 
will  be  clearly  evidenced  in  plots  to  appear  later  in  this  chapter. 


This  brief  analyses  demonstrates  the  straight-forward  manner  in  which  a 
numerical  integrator  can  be  studied  using  control  system  methodology.  The 
question  of  stability  will  be  examined  in  more  detail  later,  and  the  relation 
between  an  examination  of  the  actual  integrator  pole  locations  and  the  more 
familiar,  to  numerical  analysts,  examination  of  the  equation  for  propagation 
of  the  integrator  error  will  be  shown. 

3.2.2  The  Difference  Equation 

From  G(z)  we  can  readily  find  the  difference  equation  that  represents  the 
algorithm  to  be  programmed  on  a digital  computer.  Eq  (3-5)  yields  the  rela- 
tion 

X^(z)  = G(z)Xj(z) 

(z-l)X2(z)  = xT[yz+  (1-y)]X^(z)  (3-7) 

Since  our  objective  has  been  to  design  an  integrator,  it  is  natural  that  x^ 
be  representative  of  x and  x^  of  x.  Hence  we  say 


and  rewrite  (3-7)  as 

(z-l)X(z)  ^ \T[yz+  (1-y)]X(z)  (3-8) 

We  will,  in  fact,  designate  equality  when  we  program  Eq  (3-8).  The  measure  of 
our  efforts  will  be  how  well  x does  represent  the  true  integral  of  x.  Noting 
that  zX^  = X^^.| , we  write  the  desired  difference  equation  of  the  ZOH-TI  from 
Eq  (3-8)  in  the  following  form: 

(3-9) 


) 


n+1 


= X. 


xT[yx 


n+1 


+ (1- 
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The  ZOH-TI  is  an  implicit  integrator,  requiring  to  evaluate  • j | 

Therefore,  a predictor-corrector  procedure  or  an  extrapolation  of  x^  is  ji 

required.  The  effect  of  >■  is  to  weight  the  leading  derivative  x^^.j  versus  i 

the  lagging  derivative  x^.  Recalling  our  previous  discussion  of  the  effect 
of  V on  the  location  of  the  integrator  .dero,  we  can  intuitively  say  that 
Y ■-  '"2  tends  to  be  a lag  situation  (destabilizing)  and  > > tends  to  be  a 
lead  situation  (stabilizing).  This  lead-lag  effect  of  > will  become  very 
clear  when  we  look  at  the  frequency  response  in  Chapter  4. 

3.2.3  Phase-Shifted  Classical  Integrators  i 

Smith  has  observed  that  "many  of  the  widely  varied  classical  numerical 

integration  formulas  . . . are  actually  the  same  integrator,  differing  only 

F SI 

in  the  amount  of  phase  snift  of  the  integrand."'''  Letting  y take  on  the 

■ J 

values  0,  '2,  1 and  3/2,  we  obtain  four  familiar  integration  formulas  from 

Eq  (3-9)  as  shown  in  Table  4.  i 


Table  4;  The  ZOH-TI  and  Classical  Integrators 


Y 

Difference  Equation  (A  = 1 ) 

Name 

0 

"n-fl  = \ 

Euler 

1 

2 

Vl  = =^n  " F ^Vl  " Ni) 

Trapazoidal 

1 

1 

*n.'  “ *->  * ■'Vl 

Rectangular 

3/2 

’^n+l  " ’^n  2 ^'^_^n+l  " .^n' 

Adams  2nd  Order 

f 7 1 4. 

3.3  An  Example  Application' 

Consider  a damped  second-order  oscillator  represented  by  the  differential 
equation 

X + 2cw^x  + .O^^x  = (3-10) 


^This  is  the  same  example  as  examined  by  Smith  in  [4]  and  [S].  In  [71,  I 
explain  the  reasons  for  differences  in  our  results. 
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where  is  the  natural  frequency,  r,  is  the  damping  coefficient  and  f(t)  is  a 
forcing  function.  The  Laplace  transform  of  this  equation  is 

s^X(s)  + 2cU|^sX(s)  + tu^2x(s)  = u)^^F(s) 

We  can  draw  the  block  diagram  of  this  system  as  in  Figure  11. 


Figure  11.  Second-Order  System  Block  Diagram 


The  ZOH-TI  was  implemented  in  the  following  manner; 

1.  and  were  initially  assumed  equal  to  x^^  and  x^,  respectively. 

This  is  equivalent  to  using  an  Euler  predictor. 


2.  Eq  (3-9)  was  used  to  compute  x^^.| , which  was,  in  turn,  used  to 
evaluate  • The  process  was  iterated  until  convergence. 

An  analysis  of  block  diagrams  will  be  done  in  Chapters  5 and  6,  but  at  this 
point  it  should  be  noted  that  the  final  block  diagram  of  the  software  system 
is  not  obtained  by  simply  replacing  each  integrator  in  Figure  11  with  a ZOH-Tl. 

Tests  have  been  run  for  a number  of  forcing  functions  (including  step, 
impulse,  and  sit.usoidal)  and  the  results  are  similar  for  each.  We  will  here 
present  only  the  case  of  sinusoidal  forcing  at  resonance;  i.e.,  f(t)  = sin  i,^t. 
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this  forcing  function,  the  true  soiution  is  found  by  analytic  means  to  be 


x(t)  = !'2  e 


cos  a)  .t  + — sin  co  ,t  1 
d u) , d J 


(3-1V 


where  is  the  damped  frequency.  The  plots  that  follow  were 

computed  using  = 1 Hz,  t = 0.3,  and  >.  = 1.  The  integration  step  T and 
phase-tuning  parameter  y were  varied. 

From  Eq  (3-11)  we  note  that  the  true  solution  is  sinusoidal,  having  an 
expanding,  exponential  envelope  with  a steady  state  magnitude  of  ^ = 1.67. 
Figures  12,  13  and  14  portray  the  results  of  solving  the  differential  equation 
(3-10)  with  the  ZOH-TI  tuned  to  three  values  of  > : 0,  1,  and  ‘a,  respectively. 
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Figure  12.  ZOH-TI  Solution  for  y 
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Figure  13.  ZOH-TI  Solution  for  y = 1-0 


These  are  the  first  three  integrators  shown  in  Table  4,  and  since  it  is  well 
known  that  trapazoidal  integration  is  the  best  of  the  three,  the  results  are 
not  surprising.  The  si gni '^'i cance  of  these  data  lies  in  their  explanation 
according  to  the  laws  of  control  theory.  It  is  from  that  perspective  that  we 
gain  the  insight  as  to  what  we  have  obtained  with  the  tunable  integrator  and 
as  to  how  we  can  view  all  methods  of  numerical  integration  in  a cotmion  light. 
So  forget  for  the  moment  that  we  know  these  as  three  classical  integrators. 

All  three  integrators  perform  well  at  an  integration  step  of  T = 0.001 
sec,  and  the  curves  plotted  at  that  step  size  closely  match  the  true  solution 
of  Eg  (3-11).  At  that  step,  the  sampling  ratio  is  r = = .95  x 10"^  and 

in  Chapter  4 we  will  see  why  it  is  perfectly  reasonable  for  all  three  to 
produce  similar  results  at  that  step  size.  In  Figure  12,  > = 0,  and  we  have 


32 


o 

eri 


•O.OQ  0.50  l.DO  1.SQ  2.00  2.50  3.00  3.SQ  <(.00 


TIME 

Tigure  14.  ZOH-TI  Solution  for  y = 0.5 

intuitively  noted  that  this  situation  tends  to  produce  a lot  of  phase  lag.  In 
fact,  as  we  increase  1 the  amount  of  lag  increases  and  the  peaks  are  delayed 
and  the  solution  diverges,  having  completely  blown  up  at  T = 0.100.  (Note  the 
change  in  scale  for  that  curve.)  Quite  the  opposite  occurs  in  Figure  13  where 
Y = 1.0.  Our  intuitive  feeling  of  lead  is  buttressed  by  the  progressively 
earlier  occurrence  of  the  peaks  as  T increases.  Rather  than  going  unstable, 
however,  sensitivity  is  lost  and  the  solution  decreases  in  magnitude,  producing 
a bounded  rather  than  unbounded  error.  Note  that  the  shift  of  from  zero  to 
unity  has  completely  changed  the  sensitivity  of  the  integrator  to  increases  in 
step  size.  The  third  figure  of  the  seauence  shows  the  results  for  a balance 
between  lead  and  lag.  When  v - 0.50  vie  can  increase  the  integration  step  by 
two  orders  of  magnitude  and  still  have  a faithful  representation  of  the  true 
solution. 
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We  are  not  restricted  to  discrete  values  of  y,  however,  and  therein  lies 
the  beauty  of  tunable  integration.  In  this  example,  a of  about  0.52  appears 
to  give  the  best  fidelity.  The  continuous  shifting  of  the  solution  with 
variations  in  y is  displayed  in  Figure  15.  This  is  the  step  response  of  the 
system  and  the  true  solution  is 


x(t)  = u(t)  -e 


cos  to  .t  + ; — sin  CO  ,t 
1 d d 


(3-12) 


t = 0.5  SEC 


0.143525E1 
0.140903E1 
0.136792E1  a=.00D 
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0.133g58El 
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Figure  15.  ZOH-TI  Step-Response  Solution 


The  true  solution  at  t - 0.5  sec  is  x(0.5)  = 1.367918  and  this  is  seen  to  be 
accurately  reproduced  by  the  curve  computed  with  y = 0.5  and  T = 0.001.  The 


remaining  curves  are  computed  for  the  values  of  y Indicated  on  the  figure  ana 
T = 0.100.  For  Y=  0.53  the  solution  matches  the  true  value  to  three  decimal 
places  at  a time  step  two  orders  of  magnitude  larger  than  .001. 


i 


f 
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3.4  The  Advantages  of  Tunable  Integration 

In  this  chapter  we  have  synthesized  a discrete  approximation  to  the 
continuous  process  of  integration.  This  discrete  numerical  integrator 
possesses  the  controls  necessary  to  manipulate  its  phase  and  amplitude  charac- 
teristics in  a manner  that  allows  tuning  the  system  for  the  best  results.  The 
ability  to  tune  the  integrator  increases  its  region  of  stability  and  accuracy 
at  larger  time  steps  (see  the  discussion  in  Chapter  5 on  stability)  thereby 
allowing  for  increased  computational  speed  and  decreased  computational  cost. 
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CHAPTER  4 

Frequency-Response  Characteristics  of  the  ZOH-TI 


4.1  Why  the  Frequency  Response 

In  Chapter  2 we  saw  how  the  performance  of  a numerical  integrator  could 
be  assessed  by  comparison  of  its  magnitude  and  phase-angle  characteristics 
with  those  of  the  true  integrator  as  expressed  by  Eqs  (2-7)  and  (2-8)  and 
portrayed  in  Figure  5.  For  this  comparison  we  need  to  evaluate  the  frequency 
response  of  the  numerical  integrator.  This  has  been  done  in  [6]  for  several 
continuously  and  discretely  compensated  tunable  integrators.  The  analysis  of 
only  the  ZOH-TI  will  be  oresented  here.  Rosko^^^^  looks  at  the  frequency 
characteristics  of  a number  of  "discrete  integrators,"  such  as  Tustin  and 
Boxer-Thaler,  and  Tou^^^^  briefly  looks  at  a number  of  quadrature  formulas, 
such  as  rectangular  and  Simoson. 

There  is  a point  worth  noting  before  we  proceed,  and  it  is  implied  by 
Tou's  use  of  the  word  quadrature.  In  the  three  referenced  analyses,  the 
frequency  response  of  the  open-loop  integration  formula  is  evaluated.  The 
information  so  garnered  cannot  be  used  indiscriminantly  in  the  evaluation  of 
software  for  the  solution  of  differential  eouations,  for,  in  determining  the 
solution,  the  dynamics  of  the  differential  equation  and  the  manner  of  imple- 
menting the  algorithm  are  important.  This  open-loon  frequency-response  data 
must  be  used  with  other  information  to  get  a full  evaluation  of  the  integrator, 
but  it  is  a useful  and  important  part  of  the  overall  picture  of  integrator 
performance. 
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4.2  Methods  of  Evaluation 

There  are  a number  of  approaches  possible  to  obtain  the  desired  data,  and 
the  choice  of  which  to  use  is  relatively  free.  Three  analytic  methods  avail- 
able are  as  follows.  One  can  simply  make  the  substitution  z = e'^*^'^,  as  do 

Rosko  and  Tou,  arif^  obtain  G*(jw)  directly  from  G(z).  Smith^^^  takes  an  alter- 

1 *-w 

native  approach  to  he  same  end.  He  makes  the  transformation  z = where 
w = j tan  {uJ/2}.  (A  little  algebra  shows  that  this  relation  for  z is  equiva- 
lent to  e'^‘''^.)  The  third  analytic  approach  involves  obtaining  the  actual 
output  in  the  z-domain  bv  multiplication  of  G(z)  and  the  z-transform  of  the 
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input  and  then  taking  the  inverse  transformation  of  the  result.  Determination 
of  the  magnitude  and  phase  is  then  by  inspection;  e.g., 

A sin  (tot  + (J>)  = 1~'’  jzisin  :jt}G(z)j  (4-1) 

If  one  disdains  algebra,  the  frequency  response  can  be  determined  by 
actually  integrating  a sinusoidal  signal  and  evaluating  the  amplitude  and 
phase  of  the  resultant  output.  This  numerical  approach  is  fraught  with  a 
number  of  problems,  principally  with  respect  to  interpolation  between  data 
points  to  find  peaks  and  cross-overs  as  the  sampling  ratio  increases.  There 
is  also  a bias  generated  in  the  solution,  v/hich  must  be  subtracted  out  before 
any  evaluation  is  made.  The  third  analytic  method  is  very  tedious,  but  it 
does  enable  an  analytic  attack  on  evaluating  the  bias  in  the  numerical  compu- 
tation. The  remaining  two  approaches  are  equivalent  and  either  should  be 
preferable  to  the  novice.  The  method  used  by  Tou  and  Rosko  will  be  employed 
in  what  follows. 


4.3  Frequency  Resoonse  of  the  ZOH-TT 

To  find  G*(i(ij),  we  begin  by  writing  the  second  of  Eqs  (3-5)  in  the 
following  form: 


•,  ” i 

I -7 


Now  make  the  substitution  z 


and  obtain 


G*(j(a) 


(4-2) 


From  Euler's  equation,  we  know  that 


sin  .oT/2 


gjo)T/2  _p-ja>T/2 


Hence,  we  can  write  En  (4-2)  as 


G*(io.O 


(-jp/2) 


sin  ,T/2 
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Expanding  the  remaining  exponentials,  we  obtain  the  result 

6*(jc.)  = f (1-q)  - j cot  mT/2]  (4-3) 

This  is  written  in  a more  useful  form  by  applying  Ens  (3-6)  to  p and  g and 
noting  that 


laT  _ 0)  _ 

2 2tt/T  ""  0)^ 


r- 


The  final  result  is  then 


(4-4) 


G*(j'-o)  = Y'  1(2-, -1)  - j cot  n:|  (4-5) 

The  magnitude  and  argument  of  G*(ju>)  are  readily  found  from  Eg  (4-5) 
according  to  their  definitions  in  Egs  (2-7)  and  (2-8).  The  appropriate  rela- 
tions are 


( 2y- 1 ) ^ + cot^  r’- 1 ^ 

(4-6) 

( 2v-l  1 

(4-7) 

These  last  three  relations  tell  us  much  about  the  properties  of  the  ZOH-Tl , 

properties  which  |6]  shows  aoply  generally  to  all  of  the  tunable  integrators 

developed  to  date.  One  general  observation  is  that  the  ohase  characteristics 

of  the  ZOH-Tl  are  independent  of  X,  and  denend  on  T only  through  r.  The 

performance  of  a numerical  intearator  depends  not  on  the  freguency  of  the 

input  signal  nor  the  lennth  of  the  integration  step,  but  upon  the  relative 

size  of  that  step  compared  to  that  freguency.  A second  observation  is  that  X 

and  will  have  a scaling  effect  on  the  output  magnitude.  Thus,  we  have  a 

simple  analytic  justification  for  Smith's  observation  that  y affects  the 

f 4 1 

transient  response  and  >.  the  steady  state  solution." 

There  are  two  types  of  symmetry  apparent.  One  is  the  symmetry  of  M as  a 
function  of  , about  some  value  = 0.5.  Thus,  the  Bode  magnitude 
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Dlot  for  a rectangular  intearator  (v=l)  is  the  sane  as  for  an  Euler  integrator 
(y=0).  Smith  observes  that  this  arises  naturally  for  the  ZOH-Tl  because  "the 
zero-order  hold  integrand  reconstruction  process  introduces  exactly  a half 

1151 

period  of  lag  which  is  compensated  by  a half  period  of  lead  when  7 = + ‘2."^  ^ 
The  second  type  of  symmetry  is  the  mirror  imagery  of  as  a function  of 
r.  There  is  often,  though  not  always,  a 7^  about  which  equal  displacements  of 
Y produce  equal  displacements  of  (f.  For  the  ZOH-TI  the  value  of  y,  is  the 

V 

same  as  y^^  (i.e.,  0.5)  and  = -90°  for  all  r. 

The  last  general  property  is  the  appearance  of  the  term  cot  rn  in  all  of 
the  gain  relations.  As  r becomes  small,  this  term  dominates  the  exoression 
fo>"  and  all  integrators  have  approximately  the  same  magnitude  characteris- 
tics. If  we  expand  the  expression  for  M with  r « 1,  we  find  that  M \/^  and 
we  match  the  true  integrator  gain  slooe  of  20  db/decade.  Since  the  phase  of 
the  integrator  approaches  -90°  for  all  values  of  > when  r is  small,  we  have  a 
very  cogent  demonstration  of  why  small  integ-.-ation  steps  result  in  satisfactory 
performance  by  nearly  any  integrator. 

The  Bode  olots  of  and  i versus  log.jgr  sticwn  in  Figures  16  and  17, 
respectively,  were  generated  using  T = 0.002  sec.  Due  to  infomiation  limita- 
tions of  the  sampling  theorem,  only  values  of  r have  been  used.  The 
symmetry  of  M and  mirror  imanery  of  <p  are  evident,  as  is  the  20  db/decade 
slope  of  M for  small  r.  The  slope  remains  at  that  value  longest  for 
/ = 0.2  and  0.8,  these  then  being  the  best  choices  for  controlling  amplitude 
distortion.  It  is  interesting  to  note  that  at  r = 's  (log  r = -0.3)  the  output 
of  the  trapazoidal  integrato*^  would  be  zero  (since  the  integrand  samples  are 
being  taken  precisely  one-half  period  apart).  We  observed  in  Chapter  3 the 
difference  in  output  of  the  second  order  system  when  there  was  a lot  of  lead 
compensation  (y=1)  and  when  there  was  none  (y=0).  the  latter  case  being 
unstable  at  large  T (or  r).  Both  situations  produce  the  same  magni'-ude 
characteristics , but  orcatly  different  phase  characteristics,  for  the  open- 
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Figure  16.  ZOH-TI  Bode  Magnitude  Plot 

loop  integrator.  Since  both  gain  margin  and  phase  margin'*'  affect  system 
stability,  it  is  reasonable  to  say  that  the  phase  characteristics  are  more 
important  to  integrator  oerformance.  This  was  empirically  observed  by  Smith^ 
and  is  here  analytically  explained. 


*^Gain  and  phase  margin  are  defined  by  the  relations 

Gain  Margin  = -lG(jw^)l^j^  for  Arg  {G(ju)^)}  = -180'^ 

Phase  Margin  = 180°  + Arg  (G(juj)}  for  'G(ju),  )|  = 1 

These  terms  are  normally  applied  to  the  open-loon  transfer  function.  GH,  of  a 
closed-looD  control  system. 
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Figure  17.  ZOH-TI  Bode  Phase-Angle  Plot 


Figure  17  clearly  shows  the  lead-compensation  effect  of  > for  y > -2  and 
the  additional  lag  (beyond  -90°)  when  > < Just  as  was  observed  in  Figure 
13,  increasing  the  integration  step,  and  concomitantly  the  truncation  error, 
does  not  necessarily  induce  greater  lag  and  instability,  though  it  does  induce 
error.  The  best  tuning  to  control  phase  distortion  is  y = for  then 
'b  = -90°.  Combining  ideal  phase  characteristics  and  relatively  good  magnitude 
characteristics,  the  traoazoidal  inteorator  is  obviously,  and  well-known  to 
be,  a good  integrator.  It  is  interesting  to  note  that  except  for  y = 0.5, 
all  of  the  phase  curves  pass  through  0°  or  -180°  of  phase  at  r = 'j. 


Tuning  based  Upon  the  Frequency  Response 

The  freouency  response  data  for  the  ZOH-TI  have  been  related  after-the- 


fact  to  the  previously  discussed 


results  of  the  second-order  system  forced  at 
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resonance.  One  would,  however,  like  to  be  able  to  use  this  analysis  in  an 
a priori  manner.  For  simple  problems  the  approach  is  straightforward.  A 
tradeoff  must  be  made  between  the  competing  demands  of  amplitude  and  phase 
fidelity.  If  there  is  a clearly  dominant  frequency  in  the  integrand,  then  the 
appropriate  combination  of  sampling  ratio  and  tuning  can  be  selected  to  meet 
the  criteria  of  the  problem,  such  as  total  solution  time  or  computer  cycle 
time.  When  there  are  multiple  frequencies  of  equivalent  importance,  the  diffi- 
culty of  tuning  is  greatly  increased.  When  a single  integrator  cannot  be 
tuned  to  cover  the  frequency  band  of  interest,  then  one  possible  approach 
would  be  to  spectrally  separate  the  integrand  prior  to  integration  with 
recombination  after  integration.  The  complexities  of  such  an  approach  and  the 
difficulties  imposed  by  nonlinear  problems  have  not  been  examined  by  this 
author.  In  very  complex  problems  one  may  be  left  with  no  alternative  but  to 
tune  the  integrator  empirically.  In  any  event,  frequency  response  data  such 
as  shown  here  should  not  be  used  alone,  but  in  combination  with  other  analy- 
tical tools  such  as  the  root  locus.  This  will  be  discussed  in  the  next 
chapter. 

Return  briefly  to  our  example  problem  of  the  previous  chaoter  to  see  how 
we  might  have  proceeded  in  an  a priori  manner.  If  the  objective  is  to  maxi- 
mize the  step  size  with  phase  error  the  primary  criterion,  then  one  would 
choose  y near  0.5.  Making  y a little  larger  than  0.5  gives  a cushion  of  phase 
lead  and  also  improves  the  magnitude  characteristics.  The  values  of  r = 
for  the  four  time  steps  to  be  considered  are  0.95  x 10"^,  0.95  x 10"^, 

0.47  X 10"\  and  0.95  x 10"'  respectively  for  T = 0.001,  0.010,  0.050,  and 
0.100.  The  logarithm  of  the  largest  value  is  -1.02.  Looking  at  Figures  16 
and  17,  one  can  see  that  a value  of  y only  a little  larger  than  0.5  will 
provide  nearly  ideal  magnitude  and  phase  characteristics.  Thus,  one  should 
expect  that  y = 0.52  will  allow  use  of  T = 0.100.  The  value  of  0.52  was  used 
to  generate  the  step  response  of  Figure  18.  Compare  this  figure  with  Figure 
15,  and  note  the  increase  in  range  of  allowable  step  sizes  that  is  afforded  by 
proper  tuning.  Clearly,  freouency  response  data  can  aid  in  the  selection  of 
the  intearation  formula  to  be  used  in  the  solution  of  differential  equations. 


CHAPTER  5 

Root-Locus  and  Stability  Considerations 


r 


5.1  The  Approach 

The  ZOH-TI  has  a oole  at  z = 1 and  a zero  at  z = q = as  shown  in 

Y 

Figure  10.  Changing  the  value  of  y shifts  the  location  of  the  zero,  and  it 
also  has  a pronounced  impact  on  the  output  of  the  integrator  and  its  frequency 
response,  as  was  demonstrated  in  Chapters  3 and  4.  The  frequency  response  was 
evaluated  for  the  open-loop  integrator,  and  the  pole-zero  map  likewise  applied 
only  to  the  integrator  itself.  To  perform  a root-locus  analysis  we  must 
implement  the  ZOH-TI  in  a closed-loop  system.  How  important  the  manner  of 
implementation  is  will  be  shown  in  Chapter  6.  We  will  only  examine  a first- 
order  system  so  that  the  details  can  be  kept  to  a minimum  level  of  complexity. 
After  looking  at  this  root  locus,  we  will  discuss  the  assessment  of  stability 
using  the  root-locus  together  with  the  frequency  response. 

5.2  Root  Locus  of  a First-Order  System  with  the  ZOH-TI 
5.2.1  The  Implementation 

We  consider  the  differential  equation 

TX  + X - f(t)  (5-1) 

where  t is  the  system  time  constant.  The  block  diagram  of  this  simple  first- 
order  system  is  given  in  Figure  19. 


Figure  19.  First-Order  System  Block  Diagram’’ 


a. 

In  the  diagram  we  could  remove  the  1/-"  blocks  from  the  input  and  feedback 
paths  by  replacing  1/s  wU'b  V*s.  This  form  of  simplification  will  be  done 
in  subsequent  block  diaara  s. 

4(3 


We  have  seen  that  the  ZOH-TI  is  an  implicit  algorithm,  requiring  some 
type  of  predictor  for  or  extrapolator  for  we  shall  use  an  Euler 

predictor  in  combination  with  the  ZOH-TI  corrector.  The  section  of  computer 
code  dealing  with  the  actual  integration  of  x could  be  programmed  in  FORTRAN 
as  the  following  sequence  of  steps: 

Y = F/TAU 

DXl  = Y-X/TAU  ^ (5-2)" 

XI  = X+T*0X1  j 

DX  = Y-Xl/TAU 

X = X^-LAMDA*T*(GAMMA*DX+(l-GAMMA)*DXl) 

Equations  (5-2)  represent  the  Euler  predictor  and  Eqs  (5-3),  the  ZOH-TI 
predictor.  We  assume  that  F is  evaluated  prior  to  Eos  (5-2)  and  that  the 
appropriate  iterative  loops  are  defined. 

Let  us  look  at  these  equations  as  recursion  relations  and  attempt  to 

.J.4- 

follow  the  flow  of  information  through  the  software.  Neglect  ■f'or  the  moment 
the  association  of  any  time  base  with  the  variables,  since  we  separately  track 
the  independent  variable  on  the  digital  computer  by  incrementing  it  with  the 
integration  step.  Suppose  that  we  are  in  the  100th  pass  through  Egs  (5-2)  and 
(5-3),  having  completed  99  passes.  Subscripting  the  variables  to  indicate  the 
number  of  the  respective  evaluation,  we  can  now  write 


We  could  more  simply  write  0X1  = (F-X)/TAU,  but  have  chosen  this  form  to 
relate  to  the  subsequent  discussion  on  implementation. 

The  idea  for  this  aporoach  was  provided  by  J.M.  Smith  in  private  conversation. 
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Note  that  Xgg  is  the  stored  value  of  X until  the  last  computation. 

At  this  time,  associate  a time  base  with  X and,  recalling  that  z"^X^^.|  = X^, 
identify  the  increment  T in  time  with  a unit  increment  in  the  subscript.  Thus, 
we  can  transform  Eqs  (5-4)  into  algebraic  form: 


Y(z)  = F(z)/t 
DXl(z)  = Y(z)-z"'X(z)/-i 
Xl(z)  = z-’X(z)+T-DXl(z) 

DX(z)  = Y(z)-X1(z)/t 
X(z)  = z-'X(zU\T[yDX{z)  + (1-y)DX1(z)] 


(5-5) 


5.2.2  The  Block  Diagram 

It  is  from  Eqs  (5-5)  that  we  can  derive  the  transfer  function  of  the 
software  algorithm  and  then  perform  the  root-locus  analysis. 

The  question  at  this  juncture  is  what  is  the  appropriate  transfer  function 
for  the  integrator?  Is  it  X/DX  or  is  it  X/DXl?  Since  Xl  and  DX  are  internal 
to  the  algorithm,  we  should  seek  X/DXl.  From  the  last  of  Eqs  (5-5),  we  have 


X(z)  XT(l-v)z  XTrZ  DX(z) 

Dxit^‘" 


(5-6) 


We  can  obtain  DX/DXl  by  subtracting  the  second  and  fourth  of  Eqs  (5-5)  and 
then  substituting  for  XI  from  the  third.  Proceeding,  we  have 

DXl(z)-DX(z)  = r [z"'X(z)+T-DXl(z)-z-^X(z)]  = ^ DXl(z) 


or 


= l-V: 
DXl(z)  ' 


(5-7) 


Combining  Eqs  (5-6)  and  (5-7)  we  get  the  desired  relation 

gf-ffiy.  XT(l-vT/,)  jfy 


(5-8) 


This  same  result  could  be  obtained  from  the  methods  of  block-diagram  manipula- 
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tion,  or  signal-flow  graph  analysis,  by  drawing  the  diagram  corresponding  to 
Eq  (5-5)  and  then  reducing  it  to  simple  form. 

Equation  (5-8)  is  the  "integrator"  transfer  function.  The  integrator  we 
have,  however,  is  a result  of  the  manner  in  which  we  implemented  (or  programmed, 
if  you  will)  the  algorithm.  We  did  not  simply  substitute  Eq  (3-5)  for  the 
ZOH-TI  and  T/(z-l)  for  the  Euler  integrator  into  the  1/s  block  of  Figure  19. 

In  ^act,  we  did  not  even  program  a true  Euler  integrator  as  can  be  seen  from 
Eq  (5-5)  and  will  be  made  very  clear  in  the  next  chapter. 

In  this  formulation  we  note  that  the  position  of  the  open-loop  zero  is  no 
longer  a function  of  the  value  of  y.  The  tuning  parameters,  as  well  as  the 
ratio  T/’  (which  is  analogous  to  our  sampling  ratio  r = in  the  frequency 

response),  now  affect  the  forward  loop  gain.  As  this  gain  varies,  the  pole 
location  will  change,  thereby  tracing  out  the  root  locus. 

The  integrator  of  Eq  (5-8)  can  now  be  incorporated  in  the  overall  soft- 
ware system  by  noting  the  inputs  to  DXl  and  the  feedback  of  X given  by  Eq 
(5-5).  We  will  see  from  the  result  that  this  will  not  be  equivalent  to  simply 
substituting  Eq  (5-8)  for  1/s  in  Figure  19.  The  inouts  to  DXl(z)  are  ^ F(z) 
and  X(z)  fed  oack  through  z”Vt.  As  noted  <^or  Figure  19,  the  l/i  factor  can 
be  moved  to  the  forward-oath  side  of  the  summing  junction.  With  this  simpli- 
fication we  have  the  block  diagram  of  Figure  20. 


Figure  20.  ZOH-TI/Euler  Implementation  of  tx  + x = f 


Comparison  of  Figures  19  and  20  shows  that  the  key  difference  from  direct 
substitution  of  Eo  (5-8)  for  1/s  is  the  delay  in  the  feedback  loop  that  is 
due  to  the  Euler  predictor.  This  delay  has  a significant  impact  on  stability. 


Another  interesting  comparison  shows  the  advantage  of  using  the  ZOH-TI 
for  the  corrector  in  this  analysis.  If  we  retrace  our  steps  using  a rectan- 
gular integrator  instead  of  the  ZOH-TI,  then  the  block  diagram  of  Figure  20 
takes  the  form  shown  in  Figure  21. 


What  the  ZOH-TI  provides  is  an  ability,  other  than  by  means  of  step  size,  to 
control  the  forward  loop  gain  and  to  thereby  control  overall  system  stability. 
Note  that  in  both  diagrams  the  block  Tz/(z-l)  appears:  this  is  the  open-loop 

transfer  function  of  tfie  rectangular  integrator. 

5.2.3  The  Root  Locus 

We  initially  need  to  compute  the  ooen-  and  closed-loop  transfer  functions. 
If  the  forward  path  transfer  function  is  designated  as  G(z)  and  the  feedback 
path  as  H(z),  then  Figure  20  tells  us  that 

G(z)  = ^ (l  - £71- 

and 

H(z)  = z-' 

The  open-loop  transfer  function  is  the  product  of  these: 

GH  A G(z)H(z)  = X (1  - ^)  (5-9) 


I 
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The  complete,  closed-loop  transfer  function  X(z)/F(z)  is  given  by  the  relation 

^z|  = = rkn 

where 

K = ^ (1  - ^)  (5-11) 

Equation  (5-10)  is  the  transfer  function  of  the  computer  software  that  solves 
the  differential  equation  tx  + x = f with  an  Euler  predictor  and  a ZOH-TI 
corrector.  One  could  in  fact  very  simply  replace  the  entire  algorithm  of 
Eq  (5-2)  with  the  single  difference  equation  corresoonding  to  Eo  (5-10);  i.e., 

Vi  = 

Taking  x(t=0)  = 0,  f(t)  = u(t),  t = 1 sec,  T = 0.1  sec,  X = 1,  and  y = 0.5, 
we  find  the  values  of  x at  t = 0.1  sec  and  t = 0.2  sec  to  be  respectively 
0.095  and  0.180975  by  both  Eqs  (5-2)  and  (5-12).  The  difficulty  (or  perhaps 
advantage,  depending  on  one's  perspective)  with  Eq  (5-12)  is  the  need  to 
predict  the  forcing  function  f(t)  rather  than  the  state  x at  t^^.| . 

The  root  locus  maps  the  change  in  pole  location  as  the  forward-path  gain 
increases.  Beginning  at  the  open-loop  position  the  poles  move  toward  the 
open-loop  zeros.  Equations  (5-9)  and  (5-10)  show  a single  open  loop  pole  at 
z = 1 and  a closed-loop  zero  at  z = 0.  Thus,  the  single  pole  takes  on  the 
value  1-K,  and  as  K increases  it  moves  from  the  point  (1,0)  to  along  the 
Re(z)  axis.  When  K - 2,  the  pole  is  at  (-1,0).  For  yT/:  sufficiently  large 
we  could  make  K < 0,  but  that  would  immediately  make  the  system  unstable. 
Hence,  we  plot  the  root  locus  as  shown  in  Figure  22. 


See  ilSJ,  or  any  other  similar  text,  for  a discussion  of  sampled-data  control 
system  analysis. 

.1.1. 

Note  that  we  are  no  longer  designing  a simple  integrator,  so  that  this  pole 
at  z = 1 is  not  necessarily  'deal. 
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IM  I Z } 


"he  root  locus  can  now  be  used  tor  two  important  analyses.  The  first  is  the 
matching  of  the  pole  of  this  discrete  model  to  that  of  the  continuous  system. 
This  type  of  analysis  is  performed  in  [14].  Our  attention  in  this  report  will 
be  restricted  to  the  second  use  of  the  root  locus;  stability  analysis. 

5.3  Stability  Considerations 

5.3.1  General  Comments 

Recall  from  the  freauency-response  analysis  that  stability  does  not 
guarantee  accuracy,  and  thus  we  see  how  to  use  the  frequency-response  analysis 
along  with  root-locus  analysis.  Fidelity  of  the  integration  cm  be  specified 
by  examination  of  the  former  and  stability  of  the  result  can  be  ascertained 
from  the  latter.  One  can  also  combine  these  two  analyses  with  root-matching 
techniques,  and  then  trade  off  the  competing  reouirements. 

Before  proceeding,  note  that  stability  information  is  also  available  from 
the  Bode  plots  of  Chapter  4.  The  gain  and  phase  margins  both  being  greater 
than  zero  is  the  indicator  of  stability,  where 

Gain  Margin  = - for  Arg  {G(ju)^)l  = -180° 

(5-13) 

Phase  Margin  = 180°  + Arg  iG(ja'i)'  ^ 
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From  i^igures  16  and  17,  it  is  apparent  that  these  quantities  are  not  necessarily 
readily  found  for  the  sample  integrator.  In  these  figures,  'jj  corresoonds  to 
log  r = -3.5,  for  which  the  phase  margin  is  +90°  regardless  of  y.  The  gain 
margin,  however,  does  not  exist  for  y ^ 0.5.  For  y < 0.5  it  is  the  negative 
value  of  the  gain,  in  db,  at  r = 0.5.  These  data  alone  would  tend  to  indicate 
that  no  stability  problems  should  occur  for  any  value  of  shown.  The  integra- 
tor, however,  is  part  of  an  overall  system  and  its  lag,  contributing  to  the 
system's  lag,  can  produce  instability.  Note  also  that  T just  scales  the  gain 
in  Eq  (4-6),  thereby  only  serving  to  shift  the  magnitude  plots  along  the 
ordinate.  As  T increases,  the  curves  shift  upward  and  cause  both  gain  and 
phase  margins  to  decrease  when  y < 0.5,  thereby  decreasing  stability.  When 
y 0.5,  the  gain  margin  decreases  but  the  phase  margin  increases.  Thus,  we 
have  yet  another  bit  of  insight  as  to  what  the  impact  of  changing  step-size  is, 
insight  with  a different  perspective  than  that  of  truncation  error.  We  also 
see  why  the  frequency-response  data  should  be  used  in  conjunction  with  other 
analyses  such  as  the  root-iocus,  since  the  former  is  not  conclusive. 

5.3.2  Analysis  of  the  First-Order  System 

The  region  of  stability  in  the  z-olane  lies  within  the  unit  circle,  and 
therefore  the  samoled-data  system  represented  by  Eo  (5-10)  (which  is  our  soft- 
ware transfer  function)  is  stable  for  values  of  K in  the  interval  [0,2j.  This 
interval  presents  us  with  a range  of  possible  values  for  the  tuning  parameters 
X and  Y and  the  sampling  ratio  T/-,  all  of  which  will  yield  system  stability. 

To  see  what  choices  of  v might  be  available  to  us,  let  \ = 1;  then  the 
limits  on  K give 


0 _ K ^ 2 

0 < I ( 1 . XL  ) < 2 


(5-14) 


A plot  of  acceptable  values  y versus  T/t  then  takes  the  form  of  Figure  23. 
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Figure  23.  Stability  Region  for  Tuning 


It  is  interesting  to  note  that  the  lower  bound  on  K generates  the  upper  bound 
on  ( and  vice  versa.  This  is  quite  logical  since  larger  values  of  K move  the 
system  towards  instability  and  lower  values  of  y move  the  integrator  towards 
instability.  This  figure  also  parallels  the  results  of  Chapter  3 in  that  it 
points  out  that  for  small  values  of  T/t,  the  :hoice  of  > has  much  latitude. 
(See,  for  example.  Figures  12-14.) 

5.3.3  A Trial  Solution 

We  have  looked  at  the  frequency  response  of  the  ZOH-TI,  and  we  have  looked 
at  the  root  locus  of  the  ZOH-TI/Euler  oair  applied  to  a first  order  system  and 
the  stability  implications  of  that  aoplication.  A brief  example  demonstrates 
one  possible  approach  to  tuning  based  on  those  analyses.  Let  ^ = 1 sec, 
f = cos  t ( = 1 rad/snc)  and  choose  an  integrator  with  as  large  a time  step 

as  possible  to  give  roasonab'e  fidelity.  Figures  16  and  17  show  that  for  y in 
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the  range  .49  to  .51  we  obtain  a fairly  good  magnitooe  characteristics  up  to 
about  log  r = -1  or  r = 0.1.  Fioure  Z2  shows  this  should  be  well  within  the 
region  of  stability  of  the  system,  and  that  we  should  be  able  to  aooroach  the 
information  limit  of  T/t  = h-  l^he  value  of  r would  imply  T = — ^ = 0.63,  and 
T/"  - '-2  v/ould  imply  T = 0.5.  Let's  be  conservative  and  choose  T = 0.2,  which 
gives  r = 0.032  and  T/--  = 0.2.  From  Figure  16,  we  see  that  for  r = 0.032 
(log  r -1.5)  is  about  -40  db.  Since  that  figure  was  produced  for  T = .002 
rather  than  T = 0.2,  we  raise  the  curve  40  db  and  see  that  X should  be  about 
unity  (since  1/ia  = 1 gives  the  true  integrator  gain  as  zero  db  and  this  is 
what  we  have) . 

the  true  solution  to  Ea  (5-1)  v/hen  f(t)  = cos  ut  is 


Xf(t) 


-T_ 

1 + ( 'X  ' J 


(cos  ut  + (XT  sin  .t) 


(5-15) 


Using  X slightly  larger  than  unity,  > = 1.01,  and  y = 0.51,  the  computed 
solution  x^(nT)  compares  to  the  true  solution  as  shown  in  Table  5.  The  value 
of  the  gain  K is  0.18  f'^^om  Eg  (5-11).  Though  not  providing  excellent  accuracy, 
this  simple  combination  of  the  Euler  and  ZOH-TI  has  given  a reasonable  >'Tpro- 
duction  of  the  true  solution.  Bette>~  tuning  and  a slightly  smaller  time  step 
v/ould  improve  the  results.  Accuracy  was  not  the  point  of  this  demonstration, 
nowever.  The  objective  was  to  show  how  one  could  initiate  the  tuning  process 
by  using  tiie  frequency-response  and  root-locus  data. 


■ i. 

Note  that  one  criteria  reflects  sampling  relative  to  the  input  f(t)  and  the 
other  relative  to  the  system  Mme  constant  t.  One  should  not  exnect  two  such 
Similar  values  of  T to  arise  in  general. 

Since  T scales  the  magnitude,  as  shown  by  Eg  (^-6l,  the  curve  for  T = 0.2  is 
simply  ^‘0  db  above  that  for  T = O.DO?. 
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Table  5:  True  Solution  vs  ZOH-TI/Euler 


t 

True 

x^(t=nT) 

ZOH-TI/Euler 

Xc(nT) 

0.8 

.482 

.458 

1.0 

.507 

.473 

1.2 

.497 

.453 

3.8 

-.713 

-.705 

4.0 

-.714 

-.720 

4.2 

1 

CTi 

CO 

CO 

-.708 

^ 6.8 

.681 

.698 

7.0 

.705 

.708 

7.2 

.701 

.690 

10.0 

-.692 

-.704 

10.2 

1 

o 

-.706 

10. ^ 

-.694 

-.680 

I would  mal'e  one  more  point  before  proceeding.  I*"  the  z-transform  of 
,os  *■  is  multiplied  by  T^|^(z)  fro!?  Eq  (5-10)  and  the  result  is  expanded  by 
partial  fractions,  one  of  the  terms  will  have  the  factor  The  solu- 

p 4-  1 

tion  component  cor'^espondi ng  to  this  term  has  the  form  Cfl-K]  , where  C is 
an  explicit  function  of  K,  T and  w.  This  is  the  term  corresponding  to  the 
exponential  in  Eq  (5-15),  and  when  K > 2 it  will  oscillate  with  increasing 
size,  correspond ing  to  the  previously  identified  instability  for  values  of  K 
greater  than  two. 

5.3.4  Stability  from  Two  Perspectives 

When  most  numerical  analysts  talk  about  the  stability  of  a numerical 
procedure,  they  refer  to  the  propagation  of  errors  through  the  algorithm  and 
question  whether  or  not  an  error  introduced  at  one  point  will  die  out  or  grow 
without  bound.  Rosko  presents  a clear  presentation  of  this  type  of  analysis 
in  fl2|. 

Following  Rosko's  presentation,  take  the  multistep  method 
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(5-16) 


Ml  Mi 

* =R„ 


where  ep^  is  the  roundoff  error  at  the  step.  A similar  relation  holds  for 

the  true  value  of  x at  t , , , when  eo  is  replaced  by  er  , the  local  truncation 

n+ 1 ''n  n 

error.  The  difference  between  the  two  values  of  x is  the  error,  and  it  propa- 
gates according  to  the  relation 


m m 

3 = X)  a.e  • + T X b.e  • + E 

n+l  .‘—L  ^ n-1  ^ n-i  n 

1=0  1=1 


where  E„  = er  - en  . Given  that  x = f(x,t),  the  mean-value  theorem  allow  us 
n I n ^'n 

to  write  this  relation  as 


III 


(5-17) 


where 


Typically,  the  assumption  is  made  that  f^^  is  at  worst  slowly  varying,  so 
it  is  taken  as  a constant,  The  characteristic  equation  corresponding  to 

Eg  (5-17)  is  then,  in  terms  of  z. 


0=  .{Z)  - n-Tt..,giz"'  ' - 


(5-18) 


To  he  stable  we  want  all  the  roots  of  Eg  (5-18)  to  lie  within  the  unit  circle. 
Corisider,  for  example,  an  Euler  integrator.  The  relation  x^^.|  = + Tx^ 

gives  us  from  En  (5-16) 

m = 0 (1-step  integrator) 


^0  ^ ^ 

h = 1 
0 
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Equation  (5-18)  becomes 


0 = = z - (1  ^ Tf^)  (5-19) 

and  we  have  the  criterion  for  stability  as 

|1+TfJ<l  (5-20) 

For  real,  we  have  the  reqion  of  stability  given  by  the  relation 

Tf^e(-2,0]  (5-21) 

Compare  this  result  with  one  obtained  via  the  type  of  analysis  presented 
earlier  in  this  chapter.  For  the  first-order  system  examined,  as  we  will  see 
in  Chapter  6,  a true  Euler  integrator  can  be  implemented  to  provide  the  block 
diagram  of  Figure  24. 


Figure  24.  True  Euler  Implementation  of  tx  + x = f 


The  closed-loop  transfer  function  is 


'CL 


(z) 




z-(l-T/-) 


This  is  stable  for 


1(1-T/t)|  < 1 I 

0 < T/t  < 2 ' 


(5-22) 


(5-23) 


These  relations  are  identical  to  Eqs  (5-20)  and  (5-21),  for  f^  = -1/t. 


Now  consider  Simpson's  integrator  (see  Eq  (2-4)),  in  the  form 


n+i  n-l  3 ' n+i  n n-1 


The  parameters  of  Eq  (5-16)  are 

m = 1 (two-step  integrator) 


ao  = 0,  a^  = 1 


b 1 = 1/3,  bg  = 4/3,  b^  = 1/3 


and  the  characteristic  equation  is 

0 = = (l-T/3  f^)  - 4T/3  f^z  - (l+T/3  f^) 


(5-24) 


(5-25) 


The  roots  of  this  quadratic  are 


z 


Rosko  notes  that  for  all  real  values  of  Tf^  one  or  tne  other  of  these  two 
roots  will  be  outside  the  unit  circle,  except  for  Tf^  = 0 which  puts  the  two 
roots  on  the  unit  circle.  Thus,  as  we  oointed  out  in  Chapter  2,  th’s  inte- 
grator is  not  suitable  for  solving  differential  equations. 

If  we  look  at  the  basic  integrator  transfer  function,  as  was  done  in 
Chapter  4,  then  for  Simpson's  integrator  we  have 


G(z)  = 


73  (7^ 


z^-1 


The  poles  of  the  integrator  itself,  without  regard  to  the  system  it  is 
employed  in,  are  at  z = _^  1,  and  the  zero  outside  the  unit  circle  at  -2->’3 
will  attract  a pole  whenever  the  integrator  is  implemented  in  the  solution 
of  a differential  equation,  clearly  an  unstable  situation. 
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Thus,  in  these  two  cases  at  least,  looking  at  either  the  error  equation 
or  looking  at  the  basic  integrator  transfer  function  or  the  overall  software 
transfer  function  has  led  us  to  the  same  conclusion. 

Observe  that  for  both  the  Euler  and  the  Simpson  integrators,  the 
characteristic  polynomials  of  the  basic  inteorator  are  the  same  as  for  the 
respective  error  equation,  Eqs  (5-19)  and  (5-25),  when  the  dynamics  are 
neglected;  i .e. , 


Euler 

X 

Simoson 

^(z)!Tf._o--^=-i 

X 


The  relation  between  what  has  been  done  in  this  report  and  some  of  the  more 
familiar  approaches  to  stability  is  now  a little  clearer.  Looking  at  the 
basic  integrator  pole  locations  gives  a quick  look  at  stability  independent 
of  system  dynamics.  A feel  for  what  will  haooen  at  re‘'atively  small  time 
steps  is  gained,  for  then  Tf^  <■  1. 

The  frequency  response  data  are  based  upon  sampling  ratios  relative  to 
the  input  signal  f.  The  root  locus  analysis,  and  the  error  propagation 
analysis,  are  dependent  upon  sampling  relative  to  the  characteri sti c system 
frequencies  (recall  the  importance  of  T/t  in  our  analysis  of  this  chapter). 
Therefore,  we  have  formulated  in  these  last  three  chapters  an  initial  capa- 
bility to  design,  using  control -theory  methods,  an  integrator  to  meet 
criteria  for  sampling  relative  to  dominant  freouencies  in  both  the  driving 
function  and  the  system,  itself.  The  approach  is  not  yet  fully  or  rigorously 
defined,  but  a beginning  has  been  made. 
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CHAPTER  6 

he  Effect  of  Implementation 


6.1  Why  Worry? 

One  might  easily  be  tempted  to  take  a given  numei'ical  algorithm  and 

program  it  for  the  computer  in  any  convenient  manner  without  regard  to  the 

manner  of  programming.  That  programming,  however,  has  a significant  impact 

on  the  performance  of  the  software  and  can  make  the  difference  between  a 

[21 

successful  solution  and  an  unsuccessful  one.  Shampine,  et.  al . ' ' noted  this 
phenomena  with  respect  to  the  GEAR  and  DIFSL'E  routines  which  use  essentially 
the  same  numerical  methods. 

To  fully  analy?e  a package  as  large  as  GEAR  or  OIFSUB  in  the  manner  of 
this  chapter  would  be  a tremendous  undertaking,  and  based  upon  the  success  of 
tnose  routines  would  not  oe  justified  now.  Future  developments,  however, 
should  pay  attention,  initially,  to  the  manner  of  implementation.  The  methods 
of  control  theory  provide  a convenient  means  for  the  analysis,  for  separable 
pa»-ts  of  the  software  can  be  designed  and  their  subseouent  interactions  as 
coupled  systems  can  be  systematically  analyzed.  The  emphasis  of  this  chapter 
is  on  an  approach  to  analyzing  numerical  methods  '^ather  than  a specific 
method,  so  I will  only  present  an  elementary  examnle  of  the  considerations 
involved. 


6.?  Implenientation  of  ix  + x = f 

Let  us  take  the  same  examole  problem  as  was  used  in  Chapter  5,  but  deal 
only  with  the  predictor  part  ot  the  method  employed  therein,  the  Euler  inte- 
grator. The  way  this  in;.eqrator  was  implemented  in  Chapter  5 did  not  produce 
a true  Euler  integrator,  ’’he  difference  between  that  implementation  and  true 
Euler  implementation  is  cogently  demonstrated  by  a freouency  response  analysis 
such  as  the  one  in  Chaoter  3. 


6.2.1  A "Mixed-Mode"  Euler  Integrator 

for  convenience,  we  repeat  Eg  (5-2)  at  this  point,  but  change  DXl  and  XI 
to  DX  and  X,  resnect i vely . 


Y = F/TAU 
DX  = Y-X/TAU 
X = X+T*DX 


(6-1) 


As  noted  in  the  previous  chapter,  we  assume  that  F is  evaluated  prior  to  the 
computation  of  Y.  It  is  also  assumed  that  the  independent  variable  is  updated 
prior  to  computing  F.  Though  this  looks  like  an  Euler  integrator,  it  is  not  a 
true  Euler  integrator  because  of  the  way  in  which  DX  is  computed. 

Based  upon  the  first  oart  of  Eg  (5-4)  we  have  the  pass-count  evaluation 
(where  TI  is  the  independent  variable  time)  as  follows: 


■^T  = 

‘hoo 

TI99+T 

II 

o 

o 

u. 

0 

0 

l_ 

U- 

H 

o 

o 

>- 

F100/TAU 

'^ICO ' 

It 

o 

o 

X 

^99  "'^*^’^100 

(6-2)' 


To  compute  we  have  used  Y from  oass  100  and  X from  pass  99.  Hence, 

DX^Iqq  is  an  approximation  to  x at  some  time  which  is  correctly  associated  with 
nei ther  TIgg  nor  TI ^qq- 

In  terms  of  z,  we  convert  Eos  (6-2)  to  the  transformed  relations 

Y(z)  = F(z)/t  ^ 


DX(z)  - Y(z)  - z-‘X(z)/T 
X(z)  = z"‘X(z)  + T-DX(z) 


(6-3) 


where  we  neglect  the  first  two  relations  for  TI  and  F,  their  not  being  a 
factor  in  the  result.  From  these  relations,  we  can  write  the  integrator 
transfer  function 


■''his  is  the  actual  implementation  used  for  the  predictor  of  Chapter  5. 
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(6-4) 


This,  however,  looks  just  like  a rectangular  integrator,  having  the  corres- 
ponding ''ifferencc  equation 


n+1 


+ TDX 
n 


n+1 


(6-5) 


We  programmed  what  we  thought  was  Euler,  got  what  looks  like  rectangular,  and 
actually  have  neither  since  DX^^^  is  not  really  » combining,  os  it  does, 

^100  ^99' 

The  block  diagram  for  the  overall  system  is  similar  to  Figure  20.  The 
inputs  to  G(z)  and  the  feedback  of  X(z)  are  evident  in  Eq  (6-3).  Thus,  we 
have  Figure  25. 


Figure  25.  Mixed-Mode  Euler  Implenientation  of  ix  + x = f 


The  open  and  closed-loop  transfer  functions  are 


GH  = 


T/: 

z-1 


(6-6) 


It  is  interesting  to  compare  Figure  25  with  Figures  20  and  21,  and  to  compare 
Eqs  (6-6)  with  Eqs  (5-9)  and  (5-10).  The  forward-path  gain  is  now  simply 
K = T/  , and  the  root  locus  is  identical  to  that  shown  in  Figure  22,  with 
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the  different  definition  of  K.  Hence,  if  K = T/t  ^ 2 the  integrator  will  be 
unstable.  This  is  precisely  ^he  result  of  Eg  (5-21)  as  obtained  from  the 
error  propagation  equation.  It  is  also  the  limit  for  undistorted  reconstruc- 
tion of  a signal  imposed  by  the  sampling  theorem. 

6.2.2  A "True"  Euler  Implementation  of  tx  + x = f 
Revise  the  coding  of  Eq  (6-1)  as  follows: 

W = Y 
Y = F/TAU 
DX  = W-X/TAU 
X = X+T*DX 

What  we  have  done  is  introduce  further  delay  into  the  system  by  creation  of 
the  dummy  variable  W.  The  effect  of  the  additional  delay  can  be  seen  in  the 
develooment  of  the  closed-loop  transfer  function. 

Proceeding  as  before,  the  pass-count  relations  are 


y (6-8) 


^ (6-9) 

y 


100  " 

TI99  ^ ^ 

100  " 

>'99 

100  " 

'(^'100^ 

100  " 

Fi,g/TAU 

100  " 

*^100  ' ^99 

100  " 

X99  + T*DX 

100 


and  the  important  relations  in  terms  of  z are 


W(z)  = z*‘Y(z) 

Y(z)  = F(z)/t 
0X(z)  - W(z)  - z''X(z)/t 
X(z)  = z-’X(z)  + ^*DX(z) 


(6-7) 


‘^2 


The  integrator  transfer  function  takes  the  same  form  as  before: 

G(z)  . (6-4) 

We  do  not  have  the  same  result,  however,  because  of  what  represents.  It 

can  now  be  clearly  associated  with  because  of  the  delay  associated  with  W. 

The  inputs  to  DX  and  the  feedback  of  X provide  the  block  diagram  of 
Figure  26a,  which  can  be  alternatively  put  in  the  form  of  Figure  26b. 


Figure  26a.  True-Fuler  Implementation  of  tx  + x = f 


The  latter  figure  is  similar  to  Figure  19  and  T/(z-l)  is  the  transfer  function 
for  the  Euler  integrator.  This  is  a rare  instance  where  the  desired  integrator 
does  simply  replace  1/s  in  the  continuous-system  block  diagram. 

Open-  and  closed-looo  trans<"or  functions  for  the  true  Eule>"  implementation 

are 
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(6-10) 


Comoarison  with  Eqs  (6-6)  shows  that  the  effect  of  W has  been  to  delete  the 
closed-loop  zero  (which  had  been  a source  of  phase  lead).  The  open-  and 
dosed-looo  pole  locations  are  identical  in  the  two  fo'Tiis  of  integrator: 
z ^ (1-T/  ).  Figure  22  represents  the  root  locus  with  K = T/t.  Both  inte- 
grators go  unstable  at  T/t  = 2,  but  the  phase  characteristics  are  very  differ- 
ent orior  to  that,  as  we  will  now  see. 

6.3  Comparison  of  Frequency  Responses 

Using  the  method  of  Chapter  4 and  Eqs  (6-6),  we  find  the  sampled-data 
transfer  function,  mann^^-ude,  and  phase-angle  relations  of  the  closed-loop 
system  with  the  mixed-mode  Fuler  integrator  to  be 


Tci_(j-) 


' ■'ll-_.ll- JijJ  cos  sT  - jd-T/  i j sinmTI 
(T/r)'  +4(1-T/  ) sin'T/2 


V I 
mm 


T/-  + (1-T/O  sin-  r 


,(r)  = -tan" 


(1-T/t)  sin~  

T/>  cos  2r-  + 2 sin^  nr 


(6-11) 


Similarly,  for  the  true-Culer  implementation  we  find 


T/t  IT/t  - 2 sin'  .T/2-jsin(aT 
(T/t)'"  + 4(l-T/r)  sin''  wT/2 


( r)  = 
te ' ' 


lIZl- 


!(S  1/-V  M1-T/-)  sin'  r’T) 


t.  (ri  -tan“' 

r a ' 


sir,  ( 

T/’  - 2 sin'  r-  i 


(6-12) 
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Note  that  Eqs  (6-11)  apply  equally  to  the  ZOH-T!/Euler  combination  by  replacing 
T/-[  with  K from  Eq  (5-11).  Also  note  that  which  is  not  surprising, 

since  both  implementations  have  the  same  open-loop  transfer  function. 

Since  we  are  trying  to  solve  the  continuous-system  differential  equation, 
the  frequency  response  of  the  continuous  system  provides  a valid  criterion 
against  which  to  measure  the  two  forms  of  Euler  implementation.  The  transfer 
function  of  interest  is 


where  we  assume  (nr)^  <<  T/t.  Once  again  we  have  seen  why  at  small  integra- 
tion steps  we  have  a fairly  free  choice  of  integration  methods. 

Of  more  importance,  especially  to  the  p'"esent  discussion,  is  what  happens 
at  larger  samplina  ratios.  Consider  the  specific  value  of  T/t  = 0.3.  Then, 
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for  various  values  of  r,  we  obtain  the  data  shown  in  Table  6 and  plotted  in 
Figure  27. 


Table  6:  Frequency-Response  Comparison 


Sampl ing 
Ratio  (r) 

Continuous 
System 
Gain  (FI) 

Mi xed-Mode/ 
True-Eul er 

Continuous 
System 
Phase  ($) 

Mi xed-Mode 
Phase 

'^mm' 

T rue-Euler 
Phase 

.001 

.9998 

.9998 

- 1.200 

- 0.8399 

- 1.200 

.003 

.9980 

.9986 

- 3.595 

- 2.517 

- 3.597 

.006 

.9922 

.9945 

- 7.162 

- 5.018 

- 7.178 

.010 

.9788 

.9850 

. -11.83 

- 8.298 

- 11.90 

.030 

.8467 

.8854 

-32.14 

-22.78 

- 33.58 

.060 

.6227 

.6913 

-51.49 

-36.43 

- 58.03 

.100 

.4309 

.5018 

-64.48 

-43.49 

- 79.49 

.300 

.1572 

.2164 

-80.96 

-28.69 

-110.4 

.500 

.0951 

.1765 

-84.55 

0.0 

-180.0 

Note  that  these  data  do  not  show  what  happens  as  we  increase  T/x.  We  already 
know  that  instability  arises  for  T/x  > 2. 

■^he  difference  in  performance  to  be  expected  from  the  two  codes  with 
T/x  = 0.3  is  clear.  While  both  have  less  of  a drop  off  in  gain  beyond  the 
break  frequency  (u  - 1/:,  r = T/x/2tt)  than  does  the  continuous  system,  their 
phase  characteristics  are  dramatically  different.  The  true  Euler  implementa- 
tion has  far  more  phase  lag,  and,  if  part  of  a larger  simulation  effort,  could 
more  readily  induce  overall  system  instabil  •' ty.  The  mixed-mode  Euler  has  even 
less  phase  lag  than  the  continuous  system. 

h. 4 The  Conclusion 

The  conclusion  to  be  drawn  from  th^s  analysis  is  obvious.  One  must  not 
simply  program  an  algorithm  on  the  digital  computer  without  regard  for  the 
effects  of  the  manner  in  which  the  program  is  implemented. 
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CHAPTER  7 

Conclusions  and  Reconitnendations 


7.1  Summary 

The  fundamentals  underlying  tunable  integration  have  been  presented  in 
this  report,  and  they  have  been  analyzed  in  the  context  of  a control -theory 
approach  that  is  applicable  to  the  analysis  of  any  technique  ■''or  the  solution 
of  differential  equations. 

The  initial  discussions  were  of  a general  nature  and  defined  the  errors 
resulting  from  numerical  methods  in  terms  of  frequency- response  characteris- 
tics and  pole-zero  locations.  The  concept  of  an  ideal  integrator,  exact  for 
a specified  input,  was  used  to  motivate  the  tunable  integrator. 

A detailed  explanation  of  tunable  integv'ation  and  derivation  of  the  ZOH- 
TI  was  followed  by  a demonstration  of  tuning  for  the  solution  of  a second- 
order  differential  equation.  The  ZOH-TI  was  then  thoroughly  analyzed  for  its 
open-loop  frequency-response  properties.  The  e'ffect  of  tuning  upon  magnitude 
and  phase-angle  character! sties  was  cogently  oortrayed  in  Bode  plots  o ^"db 
and  i'.  The  use  of  these  data  for  a prior’  tuning  of  the  integrator  was 
briefly  discussed,  preparatory  to  a more  detailed  discussion  of  tuning 
following  the  root-locus  analysis. 

To  define  the  closed-loop  system  for  a root-locus  analysis,  the  first 
order  differential  equation  rx  + x = f was  used.  The  ZOH-TI  was  employed  as 
a corrector  for  a mixed-mode  Euler  predictor,  and  the  block  diagrams  and 
transfer  functions  for  tite  software  packaoe  were  developed.  These  clearly 
showed  the  effect  (and  value)  of  the  tuning  parameters  on  forward-path  gain 
and,  hence,  on  pole  position  as  the  root-locus  was  traced  out. 

A discussion  o^  stability  showed  the  coupling  between  integration  step 
size  and  the  integrator  tuning  parameters.  A stability  region  was  graphi- 
cally depicted.  The  combined  use  of  frequency-response  and  root-locus- 
stability  considerations  was  then  demonstrated  as  an  aoproach  to  a priori 
tuning.  The  final  portion  of  this  part  o^  the  report  related  the  approach  to 
stability  taken  herein  to  the  more  familiar  approach  involving  the  error- 
propagation  e((Uatior. 
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The  final  topic  dicsussed  was  the  effect  of  orogram  implementation.  It 
was  demonstrated  how  an  apoarently  innocuous  change  in  program  coding  could 
seriously  affect  the  frequency  response  of  the  integrator. 

7.2  Conclusions 

There  are  numerous  conclusions  to  be  drawn  from  the  research  reported  on 
herein.  I briefly  state  the  more  important  ones  below. 

1.  Control -theory  methods  of  analysis  are  important  and  useful  for  the 
analysis  of  numerical  methods. 

2.  The  tunable  integrator  can  be  formulated  as  a low-order  integrator 
with  an  inherent  capability  to  adjust  its  frequency-response  and  pole-zero 
characteristics.  This  provides  means  other  than  variation  of  order  or  step 
size  to  control  integrator  performance.  The  result  promises  significant 
savings  in  comuutational  burden  and  time. 

3.  Frequency-response  and  root-locus  analyses  provide  the  data  for 
initial,  albeit  elementary,  a priori  tuning  of  the  tunable  integrator. 

4.  The  method  cf  implementing  a chosen  algorithm  in  computer  code  has  a 
significant  effect  upon  software  performance.  Control  theory  provides  an 
effective  way  to  analyze  the  potential  performance  of  alternative  code  forms. 

5.  This  report  has  not.  nor  was  it  intended  to,  developed  a software 
package  that  in  any  sense  can  be  thought  of  as  a user-oriented  tool  in  the 
manner  of  packages  such  as  GEAR. 

7.3  Recommendations 

Again,  simply  enumerating,  I make  the  following  recommendations  for 
follow-on  work  in  this  a'^oa. 

1.  Proposals  have  appeared  in  the  literature  for  standardization  of  test 

probleiiis  used  in  evaluating  numerical  integrators.  Test  cases  such  as 

those  proposed  by  Krogh  and  others  employed  in  [2]  and  [19]  should  be  employed 
to  thoroughly  evaluate  the  performance  range  and  error  characteristics  of  the 
tunable  integrators. 

2.  Many  classical  integrators  can  be  shown  to  be  special  cases  of  the 
tunable  integrators.  The  theory  should  be  fullv  developed  to  tie  the  classi- 
cal and  tunable  concepts  together,  to  analyze  the  truncation  model  of  error  of 
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the  former  versus  the  phase-and-cain  model  error  of  the  latter,  and  to 
fully  exolain  the  behavior  of  tunable  integrators.  A start  toward  this  end 
has  been  made  in  the  control -theory  analysis  of  this  report. 

3.  At  present,  employing  a tunable  integrator  on  a comnlex  system 
requires  empirical  determination  of  the  proper  tuning-oarameter  values.  An 
elementary  approach  to  a priori  tuning  has  been  discussed.  Far  more  sophis- 
ticated and  reliable  approaches  are  reauired. 

4.  The  Air  Force  and  other  technological  communities  have  a potentially 
valuable  tool  at  hand.  The  capability  of  tunable  integration  to  solv  problems 
where  storage  space  and  computation  time  are  limited,  but  a degree  of  accuracy 
is  reou i red,  must  be  experimentally  determined  by  actual  implementation  when- 
ever feasible. 


1 
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